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Abstract:  This  document  presents  the  results  of  a  detailed  research  investigation  on  a 
fundamental  problem  in  wave  mechanics:  the  propagation  of  stress  waves  in  hetero¬ 
geneous  elastic  solids.  This  phenomenon  is  fundamental  to  many  disciplines  in  en¬ 
gineering  and  mathematical  physics:  seismology,  earthquake  engineering,  structure 
acoustics,  composite  materials,  and  many  other  areas.  The  theory  and  methodologies 
of  hierarchical  modeling  of  heterogeneous  materials  are  extended  to  elastodynamic 
cases  to  make  possible  the  control  of  the  modeling  error  in  the  local  average  stress. 
One-dimensional  steady  state  and  transient  applications  are  given. 
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1  Introduction 


Today,  the  use  of  composite  materials  is  wide  spread  in  industrial  and  military  ap¬ 
plications.  These  materials  generally  exhibit  a  complex  microstructure:  mechanical, 
geometrical  and  topological  features  which  are  realized  at  scales  small  in  comparison 
with  characteristic  dimensions  of  typical  structural  components.  It  is  well  known  that 
these  micro-scale  features  and  the  mechanical  properties  of  subscale  constituents  gov¬ 
ern  the  overall  response  and  service  life  of  the  structure.  Despite  this  popularity  in 
engineering  applications,  the  analysis  of  the  response  of  such  heterogeneous  materi¬ 
als  to  service  loads  is  a  very  complex  undertaking  and  has  been  the  subject  of  research 
for  several  decades. 

For  more  than  half  a  century,  work  on  the  mechanics  of  materials  has  focused  on 
determining  so-called  effective  properties:  averaged  or  smoothened  properties  that  re¬ 
flect  in  some  global  sense  the  response  of  specimens  of  the  material  to  external  loads. 
These  average  properties  are  typically  what  are  determined  by  standard  laboratory 
tests,  e.g.:  extension,  compression,  or  torsion  of  rods.  A  major  goal  of  contemporary 
research  has  been  to  determine  bounds  on  the  various  material  parameters.  Some  of 
the  most  revered  work  in  mechanics  over  several  decades  has  been  devoted  to  this 
subject.  As  well  known  examples,  the  works  of  Hill  [14],  Hashin  and  Shtrikman  [13], 
Balendran  and  Nemat-Nasser  [1],  and  Nemat-Nasser  and  Hori  [15]  are  noted.  The 
averaging  methods  also  spawned  new  mathematical  research  into  what  is  called  ho¬ 
mogenization  of  partial  differential  equations.  The  classical  works  by  Bensoussan  et 
al  [5]  and  Sanchez-Palencia  [23]  are  examples  that  provide  a  mathematical  justifica¬ 
tion  by  assuming  periodicity  of  the  microstructure.  In  more  recent  literature,  modern 
methods  of  imaging  and  computations  have  been  used  to  derive  effective  properties 
of  actual  material  specimens.  Examples  are  the  works  done  by  Ghosh  et  al  [10,  11] 
(Voronoi  Cell  Finite  Element  Method);  Fu  et  al  [8]  (Boundary  Element  Method)  or 
Terada  et  al  [25]  (using  digital  data  of  the  microstructure,  obtained  by  X-ray  Comput¬ 
erized  Tomography). 

In  a  similar  vein,  a  multi-scale  approach  has  been  used  by  Guedes  and  Kikuchi  [12] 
and  Terada  and  Kikuchi  [24],  They  use  homogenization  techniques  to  obtain  over¬ 
all  properties  for  the  macro-scale  and  include  asymptotic  corrections  to  account  for 
micro-mechanical  effects  (periodicity  of  the  microstructure  is  assumed).  Ghosh,  Lee 
and  Moorthy  [9]  use  a  similar  multi-scale  technique  to  study  elasto-plastic  material 
behavior.  Fish  et  al  [6,  7]  apply  a  multi-scale  technique  to  problems  of  wave  prop¬ 
agation  through  heterogeneous  materials  as  well,  but  enhance  the  method  by  also 
introducing  multiple  temporal  scales  to  capture  the  dispersion  phenomenon  caused 
by  the  material  inhomogeneity. 

All  the  previous  approaches  can  be  characterized  by  their  restriction  to  materials 
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with  periodic  micro-structures.  In  recent  years,  a  completely  new  line  of  research  on 
heterogeneous  materials  has  emerged  in  which  full  account  of  micro-mechanical  fea¬ 
tures  of  materials  can  be  made  in  predicting  macro-mechanical  behavior.  This  area  is 
referred  to  as  Hierarchical  Modeling  and  involves  the  use  of  only  enough  micro-scale  in¬ 
formation  to  determine  essential  features  of  the  macro  response  to  within  preset  levels 
of  accuracy.  Homogenization  is  used  only  as  a  step  in  a  broader  algorithm.  Oden  and 
Zohdi  [21]  and  Zohdi  and  coworkers  [29,  28]  introduced  a  hierarchical  adaptive  mod¬ 
eling  method  for  problems  in  elastostatics  based  on  global  error  bounds.  The  method 
is  aimed  at  providing  a  hierarchy  of  descriptions  of  the  physics  (or  scales),  that  can 
be  used  in  different  subdomains  of  the  material.  Instead  of  heuristically  choosing  a 
level  of  description  for  each  subdomain,  a  mathematical  tool  in  the  form  of  a  posteriori 
error  estimates  of  the  modeling  error  is  used  to  identify  what  level  of  sophistication  is 
needed  in  each  subdomain. 

The  original  method  was  referred  to  as  the  Homogenized  Dirichlet  Projection  Method 
(HDPM).  It  involves  two  levels  of  descriptions:  a  homogenized  macro  description 
and  the  exact  micro-mechanical  description.  The  algorithm  proceeds  as  follows:  ini¬ 
tially  homogenized  material  properties  are  used  throughout  the  entire  domain;  next, 
a  posteriori  error  estimates  of  the  modeling  error  are  obtained  and  an  iteration  process 
is  started  in  which  in  critical  regions  of  the  material  the  fine-scale  problem  is  solved 
by  using  the  homogenized  solution  as  Dirichlet  boundary  data  on  the  boundary  of 
the  subdomain;  the  iteration  process  continues  by  including  more  and  more  critical 
regions  into  the  fine-scale  analysis  until  the  error  estimate  meets  certain  user-preset 
tolerances. 

Where  the  previous  procedure  uses  global  error  estimates,  Oden  and  Vemaganti 
[19,  20,  26,  27]  advanced  this  work  by  introducing  the  Goal-  Oriented  Adaptive  Local 
Solution  Algorithm  (GOALS),  where  error  estimates  in  local  quantities  of  interest  are 
used.  Such  estimates  allow  one  to  use  goal  oriented  adaptive  strategies,  in  which  a 
model  is  adapted  so  as  to  yield  accurate  characterization  of  specific  features  of  the 
response  identified  by  the  analyst. 

As  mentioned  before,  these  works  were  all  done  within  the  framework  of  elas¬ 
tostatics.  Recently,  an  extension  of  the  GOALS  philosophy  to  general  goal  oriented 
engineering  applications  has  been  proposed  by  Oden  and  Prudhomme  [17, 18].  Their 
approach  entails  a  residual-based  analysis  of  the  modeling  error  and  has  been  inspired 
by  the  work  of  Becker  and  Rannacher  [2]  on  goal  oriented  estimation  of  discretization 
errors. 

This  investigation  presents  an  extension  of  the  GOALS  philosophy  to  the  elasto- 
dynamic  problem  by  extending  the  existing  theory  to  cover  complex-valued  solutions 
and  sesquilinear  forms  encountered  in  frequency-domain  formulations  of  the  wave 
equation.  The  notations  and  the  model  problem  are  first  introduced  in  Chapter  2.  A 
detailed  analysis  of  the  modeling  error  for  the  wave  equation  in  the  frequency  do¬ 
main  is  then  given  in  Chapter  3,  which  is  subsequently  incorporated  in  Chapter  4. 
There,  an  adaptive  modeling  technique  is  proposed,  representing  the  extension  of  the 
GOALS  algorithm  to  the  elastodynamic  case.  The  latter  two  chapters  also  include  one¬ 
dimensional  steady  state  and  transient  numerical  applications  and  verifications  of  the 
new  methodology.  Lastly,  Chapter  5  lists  some  concluding  remarks. 

It  is  noted  that  part  of  this  work  is  also  contained  in  [22], 
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2  The  Elastic  Wave  Problem 


In  Section  2.1,  the  notations  and  the  model  problem  of  elastic  wave  propagation  are  in¬ 
troduced.  Subsequently,  the  weak  formulation  of  the  wave  equation  in  the  frequency 
domain  is  defined  in  Section  2.2. 


2.1  Model  Problem  and  Notations 

Let  SlcR2  be  a  bounded  open  domain  with  Lipschitz  boundary  dQ,  and: 

dn  =  ruurt,  runrt  =  0, 

where  F„  denotes  the  part  of  the  boundary  with  prescribed  displacements  U,  and  F  t 
the  part  subjected  to  traction  and  damping  conditions.  The  domain  S2  is  the  interior  of 
an  elastic  solid  material  with  a  microstructure  that  exhibits  highly  oscillatory  material 
properties  (see  Figure  2.1).  Since  we  assume  each  of  the  constituents  to  be  linearly 
elastic,  the  relation  between  the  stress  tensor  a  and  strain  tensor  e  is  governed  by  the 
following  linear  relation: 


ct  =  Ee,  (2.1) 

where  E  =  E(x)  €  Loc(0);V2x,v^  denotes  the  fourth  order  elasticity  tensor,  which  satis¬ 
fies  the  following  symmetry  and  ellipticity  conditions: 


Eijki  (x)  —  Ejiki(x)  —  Eijik(x)  —  Ekiij  (x) , 

£,ij£,ij  —  ^ijkl{^)Cij^,kl  —  £ij£iji 


1  <  i,j,  k,  l  <  3,  ao,a\  G  M,  no  >  0,  a±  >  0, 

£ij  =  £ji  G  M2x2,  for  x  E  R2,  a.e. 

The  deformations  in  the  material  are  assumed  to  remain  small,  such  that  the  strain- 
displacement  relationships  are  linear: 

e  =  I(Vu  +  (Vu)T).  (2.2) 

In  this  analysis,  stress  wave  propagation  is  considered  only.  This  means  that  both  the 
stresses  and  displacements  are  continuous  in  fh  Now,  at  t  =  0  the  initial  displacement 
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Figure  2.1:  The  model  problem. 


and  velocity  fields  are  given  by  Uo  G  (f/1(n))2  and  Vo  G  (H1^))2,  respectively  The 
linear  elastodynamic  model  problem  can  then  be  formulated  as  the  following  initial 
boundary  value  problem  (IBVP): 


For  x  G  n,  t  >  0,  find  u(x,  t)  such  that: 


P(x)u(x,f)  —  V  •  (E(x)Vu(x,t)) 

=  f  (M), 

in  U, 

EVu(t)  •  n  +  /3u(t) 

=  t(f), 

on  Tt  Vt  >  0, 

u(t) 

=  U(f), 

on  Tm  Vt  >  0, 

u(x,  0) 

=  U0(x), 

in  Q  at  t  =  0, 

u(x,0) 

=  V0(x), 

in  FI  at  t  =  0, 

(2.3) 


where  the  function  f  G  (L2(Jl))2  represents  the  body  forces,  p  G  p>  0  a.e.  in 

the  mass  density  distribution  of  the  material,  and  f3  G  L°°(Tt)  the  damping  coeffi¬ 
cient.  Notice  that  the  tractions  t  are  in  L2(Tt)  and  that  Uo  and  Vo  are  assumed  to  be 
consistent  with  the  boundary  data. 

In  Chapter  4,  an  adaptive  modeling  technique  is  proposed  that  solves  the  wave 
problem  in  the  frequency  domain.  Toward  this  purpose,  we  apply  a  classical  Fourier 
transformation,  defined  as  follows: 

/•OO 

u(x,ca)  =  JT(u(x,f))  =  /  u(x,t)  e~,wt  dt, 

Jo 

where,  i2  =  —  1,  to  represents  the  radial  frequency,  and  u(x.  uj),  the  Fourier  transform 
of  u(x,  t),  is  a  complex-valued  function.  The  equivalent  formulation  of  (2.3)  can  then 
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be  recast  as: 


For  every  u>,  find  u(x,  uj)  such  that: 


—pu2u(x,Lu)  —  V  •  EVu(x,w) 


f(x,ta)  +  i<jjp\Jo(x ) 
+pV0(x),  in  If, 


EVu(w)  •  n  +  / 3uiu{u ) 
u(w) 


t(cj)  +  /3U0,  on 
U(w),  on  r„. 


(2.4) 


where  f,  t,  and  U  denote  the  Fourier  transforms  of  f,  t,  and  U,  respectively  Upon  solv¬ 
ing  (2.4)  for  every  frequency,  the  solution  u(x,  t)  is  retrieved  by  applying  the  inverse 
Fourier  transformation: 

1  t'°° 

u(x,  t)  =  JF-I(u(x,u;))  =  —  /  u(x,cj)  elut  du>. 

2tt 


2.2  The  Weak  Formulation 

The  space  of  complex-valued  test  functions  V  is  defined  as  follows: 

V={ve(H1(tt))2:  v,ru  =  0}  .  (2.5) 

The  equivalent  variational  formulation  of  the  model  problem  (2.4)  is  governed  by: 


For  every  u  find  u  €  {U}  +  V  such  that : 

£>(u,  v)  =  £(v),  Vv  €  V, 

where  the  sesquilinear  form  £>(., .)  and  linear  form  £(.)  are  defined  as  follows: 

C  :  V — >C,  B  :  VxV — ►  C, 


(2.6) 


£(v) 

=  /{f  +  ^U0  +  pV0}:vdx+  J  (t-F/fUo)  :  vds 

B(  u,v) 

=  «4(u,  v)  —  lo2C(vl,  v)  +  iuV(  u,  v), 

„4(u,v)  =  /  EVu  :  Vvdx, 

Jn 

(2.7) 

C(u,v)  =  /  pu  :  vdx, 

Jn 

P(u,v)  =  f  /3u  :  vdx, 

J  rt 

where  v  denotes  the  complex  conjugate  of  v.  The  sesquilinear  forms  A(-,  •),  C(-,  •),  and 
£>(•,•)  incorporate  respectively  the  elastic  deformation,  inertia,  and  damping  condi¬ 
tions  of  the  elastic  body  into  the  variational  formulation. 
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3  Modeling  Error  Estimation 


In  this  chapter,  estimates  are  derived  of  the  modeling  error  in  the  local  average  stress,  in¬ 
duced  by  the  employment  of  approximate  material  models  for  the  elastic  wave  equa¬ 
tion  in  the  frequency  domain.  These  estimates  are  employed  in  the  adaptive  modeling 
algorithm  in  Chapter  4  to  assess  the  modeling  error. 

Section  3.1  begins  with  definitions  of  the  exact  and  approximate  variational  for¬ 
mulations  of  the  wave  problem  in  the  frequency  domain.  Subsequently,  in  Section  3.2, 
the  variational  problems  are  presented  that  govern  the  corresponding  modeling  er¬ 
rors.  Elliptic  representations  of  the  error  functions  are  derived  in  Section  3.3,  which 
are  first  used  in  Section  3.4  to  derive  upper  and  lower  error  bounds.  Secondly,  a  pos¬ 
teriori  global  and  local  error  estimates  are  derived  in  Section  3.5.  One-dimensional 
numerical  verifications  of  the  error  estimates  and  bounds  are  shown  in  Section  3.6. 


3.1  Exact  and  Approximate  Material  Models 


The  quantity  of  interest  is  defined  in  functional  form  and,  since  the  wave  problem 
is  solved  in  the  complex  frequency  domain,  its  equivalent  representation  is  needed 
in  this  setting.  Thus,  if  Q(u)  £  M  denotes  the  quantity  of  interest  in  the  real-valued 
time  domain,  then  its  equivalent  in  the  complex- valued  frequency  domain,  Q( u),  is 
governed  by  the  following  transformations: 


W 


T 


V 


Q 


Q 


T 


C 


Q  =  To  (Qo T~ 1 ) , 


where  W  denotes  the  space  of  admissible  displacements  in  time  domain.  In  this  work, 
the  quantity  of  interest  is  the  average  stress  over  a  small  subdomain  S  C  If  in  the  direc¬ 
tion  of  a  user-specified  vector  w. 

Q(u)  =  ]^|  ^(Wu  •  ^) dx-  (3.1) 

Thus,  the  primal  and  dual  problem  [17]  for  the  exact  description  of  the  material  model, 
are  given  by: 


B(fl,v)  =£(v), 

Vv  £  V, 

primal  problem 

£(w,p)  =  Q(w), 

Vw  £  V, 

dual  problem 

(3.2) 
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where  p  denotes  the  influence  function  [17]  with  respect  to  the  quantity  of  interest 
Q(-),  and  the  functionals  £>(•,  •)  and  £(•)  are  defined  in  (2.7). 

Now,  the  approximate  problem  is  posed  by  using  an  approximate  description  of 
the  material  model,  characterized  by  the  functions  Eo(x)  and  po(x)-  These  functions 
represent  approximations  of  the  actual  distribution  of  the  modulus  of  elasticity,  E(x), 
and  mass  density,  p(x). 

The  approximations  Eo(x)  and  po(x)  may  represent  averaged  or  homogenized  fields 
and  are  expected  to  be,  in  general,  smoother  than  E(x)  and  p(x).  In  later  sections,  more 
explicit  definitions  of  these  functions  are  given.  At  this  point,  Eo(x)  and  po(x)  are  used 
to  define  the  following  functionals,  analogous  to  those  in  (2.7): 

#o( U0,v)  =  -40(uo,v)  -  w2c0(u0,v)  +iujV(u0,w), 

A0(u0,v)  =  J  E0Vu0:Vvdx,  (3.3) 

C0(u0,v)  =  /  p0u0  :  vdx. 

Jn 

The  approximation  of  (3.2),  can  now  be  introduced  as: 


(3.4) 


Note  that  in  the  definition  of  the  functional  the  exact  representation  of  the 

damping  condition  (the  functional  V(-,  •))  is  incorporated.  In  subsequent  applications, 
this  condition  acts  only  on  a  small  part  of  the  boundary  and  can  be  easily  implemented 
with  low  computational  cost. 


3.2  The  Modeling  Errors 

The  (modeling)  error  functions  §0  and  £0  for  the  primal  and  dual  problem,  respec¬ 
tively,  are  defined  as  follows: 

e0  =  u  —  u0,  £0  =  p  — p0.  (3.5) 

Furthermore,  the  residual  functionals  characterizing  the  accuracy  of  the  approximate 
solutions,  assume  the  forms: 


ft(uo,v)  =  £(v)  -£(u0,v),  veh, 
7£(p0,w)  =  Q(w) -6(w,p0),  weF 


(3.6) 


Now,  by  substituting  the  variational  approximate  problems  of  (3.4)  into  the  two  above 
expressions,  explicit  expressions  for  the  residual  functionals  in  terms  of  the  approxi- 
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mate  solution  pair  (uo,po)  are  obtained: 


•  Vv  —  ptu2jo  uq  •  v}  dx, 


{ETqVw  •  Vpo  -  puj2jo  W  •  p0}  dx. 


(3.7) 


Here  the  deviation  tensor  1$  and  function  jo  are  defined  in  the  following  manner: 

J0  =  I-E”1E0,  j0(x)  =  l-^.  (3.8) 

PW 

Conversely,  by  substituting  the  exact  formulation  (3.2)  into  (3.6)  and  by  employing  the 
sesquilinear  property  of  £>(•,  •),  the  following  set  of  variational  problems  that  govern 
the  error  functions  §o  and  Sq  are  derived: 


J33 

ro> 

o 

II 

^o,v), 

Vve  V, 

#(w,e0)  = 

^(po,w), 

Vw  G  V. 

(3.9) 


By  applying  Theorem  1  in  the  work  by  Oden  and  Prudhomme  [17],  and  noting  that 
both  Q(-)  and  £>(■,•)  are  linear  functionals,  the  following  expression  is  obtained  for  the 
error  in  the  quantity  of  interest: 

Q(e0)  =  72(u0,p0)  +  ^[ft(u0,£0)  +  ft(p0,e0)]. 

Substitution  of  (3.9)  finally  yields: 


Q(e0)  =  7e(u0,po)+S(e0,£o)- 


computable 


(3.10) 


In  general,  the  first  term  in  the  RHS  is  a  computable  term,  whereas  the  second  is  not. 
In  the  case  of  elastostatics,  Oden  and  Vemaganti  [19,  20,  26,  27]  accurately  estimate 
the  unknown  term  by  using  an  advantageous  property  of  the  corresponding  bilinear 
functional  in  that  case:  it  defines  an  inner  product  on  the  space  of  admissible  test  func¬ 
tions.  This  enables  the  estimation  of  the  unknown  term  by  global  norms  of  the  error 
functions,  which  subsequently  can  be  estimated  in  terms  of  the  computed  approxi¬ 
mate  solutions.  However,  due  to  the  minus  sign  in  front  of  C(-,  •)  and  the  presence  of 
the  damping  term  T>(-,  •)  in  the  definition  of  £>(•,  •),  this  sesquilinear  form  does  not  de¬ 
fine  an  inner  product  on  V  x  V  and  the  approach  proposed  by  Oden  and  Vemaganti 
cannot  directly  be  applied  to  the  elastodynamic  problem. 

To  overcome  this  difficulty,  elliptic  equivalents  of  the  error  functions  are  intro¬ 
duced  in  Section  3.3,  which  serve  as  an  intermediate  step  in  estimating  the  unknown 
term  in  (3.10). 
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3.3  Elliptic  Error  Representations 


Elliptic  representations  e  i  and  £\  of  the  primal  and  dual  error  functions  e()  and  £q  are 
defined  as  follows: 


W(Si,v) 

> 

o 

<01 

II 

=  K{  u0,v), 

Vv  e  V, 

TL{ w,ei) 

=  B{w,£0) 

II 

SI 

O 

Vw  €  V, 

(3.11) 


where  the  sesquilinear  functional  H is  defined  by: 


H:  VxV — >C, 

H(  u,  v)  =  *4(u,  v)  +  lo2C(  u,  v), 


(3.12) 


where  the  functionals  A( •,  •)  and  C(-,  •)  are  given  in  (2.7).  It  is  observed  that  the  above 
functional  differs  with  £>(•,  •)  in  the  plus  sign  in  front  of  the  inertial  term  C(-,  •)  and  the 
elimination  of  the  damping  term  £>(•,•)• 


Remark  3.3.1  The  sesquilinear  form  H(-,  •)  is  positive  definite,  hermitian,  and  consequently 
defines  an  inner  product  on  V  x  V.  The  norm  that  is  implicitly  defined  through  the  inner 
product,  is  then: 


r— - r  HV,W) 

IM|h  =  \/H(v,v)  =  sup  — n — T - .  (3.13) 

wev  I!wI!h 

By  recalling  expression  (3.10),  the  error  can  be  rewritten  in  terms  of  the  elliptic  repre¬ 
sentations  by  applying  the  variational  formulation  (3.11)2,  which  yields: 

Q(eo)  =  72(u0,po)  +  H(e0,£i). 

By  combining  (2.7)  and  (3.12),  the  last  term  in  the  RHS  can  be  expanded  as  follows: 

H(e0,£i)  =  £(e0,  <§i)  +  2w2C(e0,  £i)  -  if). 


Thus, 


Q(e0)  =  72(u0,po)  +  S(e0,£i)  +  2w2C(e0,ei)  -iujV(e0,i i). 

Substitution  of  (3.11) 1  then  gives: 

Q(e0)  =  n(u0ip0)  +  n(e1,£i)  +  2io2C(e0,£i)-icoV(e0,ii).  (3-14) 

To  estimate  the  error  Q(§o),  as  given  in  (3.14),  the  following  intermediate  error  esti¬ 
mator  7*  is  proposed: 


7*  =  72(u0,p0)  +  H(ei,£i), 


(3.15) 


where  the  terms  in  (3.14)  have  been  neglected  that  involve  eo- 
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Remark  3.3.2  At  this  stage,  no  theoretical  proof  of  the  accuracy  of  7*  has  been  given,  but 
numerical  verifications  in  Section  3.6  show  that  this  estimator  provides  a  remarkably  accurate 
estimate  of  the  modeling  error  in  the  local  average  stress,  evert  for  high  frequencies  and  for  a 
wide  range  of  material  properties. 

In  practical  applications,  of  course  one  does  not  want  to  use  7*  as  an  estimate  of 
the  modeling  error,  as  it  involves  the  solution  of  (§1,  ei).  These  functions  are  governed 
by  the  variational  problems  in  (3.11)  and  involve  the  exact  description  of  the  material 
model.  Here,  the  estimator  7*  is  introduced  as  an  intermediate  step  toward  deriving 
computable  estimators.  Numerical  experiments  show  that  7*  represents  an  accurate 
error  estimate.  Thus,  accurately  estimating  7*,  indirectly  leads  to  an  accurate  estimate 
of  the  modeling  error  itself. 

The  inner  product  property  of  Ti (•,  •)  can  now  be  used  (see  Remark  3.3.1).  Accord¬ 
ing  the  Polarization  Formula  [16],  the  last  term  in  (3.15)  can  be  expanded  as  follows: 


W(ei,ei)  =  ^||ei  +£i\\n  -  ^||Si  -  £i||h  +  |||Si  +  iii\\n  -  |||Si  -  ii i\\n- 


(3.16) 


Hence,  (3.15)  can  be  rewritten  as: 


7*  =  ^(u0,p0)  +  ^l|ei  +  ii\\n  -  ^||Si  -  £i|| n 


(3.17) 


Now,  since  the  first  part  in  the  RHS  is  computable,  estimates  and  bounds  on  7*  can 
be  derived  by  estimating  and  bounding  the  remaining  terms  involving  the  norms  of 
§1  and  i\.  In  the  following  two  sections,  a  selection  of  bounds  and  estimates  of  these 
norms  is  derived  and  used  to  propose  estimators  and  bounds  of  7*  and,  therefore, 
indirectly  of  the  modeling  error  in  the  average  stress  Q(e 0)  itself. 


3.4  Computable  Upper  and  Lower  Error  Bounds 

In  this  section,  upper  and  lower  bounds  on  the  error  estimator  7*  are  derived.  As  men¬ 
tioned  in  the  previous  section,  it  is  assumed  that  these  bounds  hold  for  the  modeling 
error  itself  as  well. 

Lemma  3.4.1  Let  (uo,po)  and  (§i,£i)  be  solutions  to  (3.4)  and  (3.11),  respectively,  and  let 
the  deviation  tensor  and  function  To  and  jo  of  (3.8)  be  piecewise  continuous  functions  in  Li. 
Then  the  following  computable  upper  and  lower  bounds  hold  for  the  norms  in  (3.17): 

||Si  +  ei||H<»7+p’ 

hiow  —  II®1  ~~  hIIw  —  Vuppi 

(3.18) 

r^<  \\ei  +  iii\\n<riupp,v 

7lloW,i  —  II®1  _  ^illw  —  huppp ! 
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ivhere: 


hupp  =  \J l|2o(uo  ±  Po)  Il3i  +  w2||  jo  (uo  ±  Po)  He, 


V 


± 

upp,i 


sj |jJo(uo  ±  *Po)  11^  +  w2||io  (Uo  ± l, 


,  ±  =  |^(uo±po,u0  +  6i::i=po)| 

Vlow~  K  +  ^po||h 

±  =  |fc(uo±fpo,uo  +  fl.fpo)| 

ll0^~  j|u0  +  Ofpo\\n 


(3.19) 


l^  =  ^(v,v) 


v||g=C(v,v) 


ivhere  the  sesquilinear  forms  and  C(-,-)  are  defined  in  (2.7),  the  residual  functional 

12(-,  •)  is  defined  in  (3 .7)1,  and  where  9±  and  Of  €  C. 

Proof:  Only  the  proof  for  the  bounds  on  the  first  norm  in  (3.18)  are  presented.  The 
proofs  for  the  other  norms  follow  analogously.  Applying  the  definition  of  the  norm 
|| -||ft  as  given  in  (3.13),  gives: 


ei  +£i||w  =  sup 
wev 


\njep  +  £i,w)| 

IHIw 


By  recalling  Remark  3.3.1  that  H(-,  •)  is  sesquilinear  and  hermitian,  this  expression  is 
rewritten  as: 


ei  +ii\\n  =  sup 
wen 


|H(ei,w)  +?f(w,£i)| 

IHIw 


and,  subsequently,  the  variational  problems  given  in  (3.11)  are  introduced,  which 
yields: 


§1  +£i\\h  =  SUP 

wev 


H(uq,w)  +  72(p0,w)| 

IHIw 


If  the  explicit  expressions  for  the  primal  and  dual  residual  functional  in  (3.7)  are  com¬ 
pared,  it  is  observed  that  7£(po,  w)  =  1Z(po,  w).  Hence, 


ei  +  £i|lw  =  sup 
wen 


|7e(u0  +p0,w)| 

IHIw 


(3.20) 


The  lower  bound  r/^w  now  follows  by  applying  the  definition  of  the  supremum  and 
choosing  w  =  uq  +  0  1  p0,  where  O'  G  C  is  arbitrary.  Thus, 


ei  +£i|| w  > 


|72(uo  +  po,u0  +  6'+po)| 

||u0  +  6>+p0||w 


V0+  eC. 
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To  prove  the  upper  bound,  (3.20)  is  rewritten  by  using  (3.7) 1  and  (2.7),  and  by  recalling 
that  Tq  and  jo  are  piecewise  continuous,  which  gives: 


|ei  +£i||  h  =  sup 
wev 


|*4(X0(u0  +  po),  w)  -  uj2C{j0{u0  +  po),  w)| 


I|W||  n 


(3.21) 


Applying  the  triangle  and  Schwarz  inequalities  to  the  nominator  of  (3.21),  gives: 

|^(T0(u0  +  po),w)|  +w2|C(j0(u0  +  po),w)|  < 


V  "d(Jo(uo  +  Po)  ,^o(uo  +  Po))  V  -d(w’ 


w 


+u2  sfc  ( jo  (u0  +  Po ) ,  jo  (u0  +  Po ) )  VC  (' w  > 


w 


<  \[m?o^o  +  Po),2To(uo  +  Po))  +  w2C(jo(u0  +  po)Jo(uo  +  Po)) 

x  yM(w,w)  +  ca2C(w,w) 

=  \J li^o(uo  ±  Po)H3v  +  w2llio  (u0  ±  Po) lie  llwll«- 

Finally,  the  upper  bound  ?/+  is  obtained  by  substituting  this  expression  into  (3.21).  ■ 


To  obtain  near  optimal  values  of  the  complex  numbers  9±  and  Of,  a  procedure  sum¬ 
marized  in  Appendix  A  is  employed.  If  the  upper  and  lower  bounds  of  this  lemma 
are  introduced  to  (3.17),  upper  and  lower  bounds  on  the  real  and  imaginary  parts  of 
7*  follow  automatically. 

Corollary  3.4.1  Given  the  bounds  riupp(  %)  am ^  Tlk>w(  i)>  ^ e  rea l  am d  imaginary  parts  of  7*  are 
bounded  asfollozvs: 


hlow,real  —  {7  }  hupp, real  1 


Vlmv,imag  XX?7T,{7  }  X  hupppmagi 


(3.22) 


inhere: 

Viow.real  =  K e  {^(u0 ,  po) }  +  -  \rf^ VP2, 

1  2  1  2 

hupp, real  =  Tie {ft(u0, p0)}  +  3r/+p  ~  jrjfow  , 

(3.23) 

Vlow,imag  =  Irn {72(u0, Po)}  +  I Vum*  -  I Vupp/> 

VuppMag  =  Zrn {72(u0, p0)}  +  \ntpp,i  ~  bfow,i  • 

Proof:  Only  the  upper  bound  on  Tie  {7*  }  is  proved,  as  the  proofs  for  the  other  bounds 
follow  analogously.  Taking  the  real  part  of  (3.17),  gives: 

■Re{y*}  =  ne{n{uo,po)}  +  j\\ei+ii\\h-^\\ei~£i\\h. 
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By  applying  the  inequality  of  (3.1 8) 1 ,  one  obtains: 

ne{^*}  <  TZe{n{uo,po)}  +  \vuVV2  -  \\\*i  -  ii\\n- 

Finally,  by  substituting  (3.18)2,  the  assertion  is  established: 

^e{7*}  <  Ke{K{ u0,p0)}  +  \vupp2  ~  • 


3.5  Modeling  Error  Estimates 

In  this  section,  two  types  of  a  posteriori  estimates  of  the  modeling  error  are  derived. 
Section  3.5.1  introduces  estimates  of  the  modeling  in  the  local  average  stress  in  a  sub- 
domain  S  of  O,  whereas  Section  3.5.2  presents  global  error  estimates  in  the  norm  j  •  1 1  . 


3.5.1  Local  Error  Estimates 


Numerical  experiments  confirm  that  the  bounds  on  the  modeling  error  of  (3.22)  in¬ 
deed  are  bounds  to  both  7*  and  Q(§o).  However,  for  the  elastostatic  case  [19],  these 
types  of  bounds  are  a  factor  of  100  higher  and  lower  than  the  error.  On  the  basis  of 
extensive  numerical  experiments  it  has  been  observed  that  for  the  elastodynamic  case, 
the  bounds  of  (3.22)  are  a  factor  1000  or  10,000  higher  and  lower  than  the  error. 

These  bounds  are  distributed  roughly  equally  around  the  error.  Accordingly,  the 
first  estimator  7avg  that  is  proposed,  uses  the  averages  of  these  bounds.  Thus,  given 
the  bounds  7low,real/  ^/low,imag/  Nipp.rea l >  <irid  r/upppmag  of  (3.22)  and  (3.23),  an  estimator 
7avg  of  7*  is  introduced  as: 


Tavg  2  |  (^/upp.real  T  ??low,real)  T  f(7upp,imag  T  l/low.imag)  J' 


(3.24) 


Numerical  tests  (see  Section  3.6),  also  suggest  that  the  bounds  7/Lypp,  rj^  and  r/^w, 
r/P)w  .  have  good  effectivity  indices  with  respect  to  the  norms  |e  i  ±  £i\\n  and  ||  ep  ±  i£\  \  \n 
they  bound.  The  upper  bounds  7^  and  rj^  i  exhibit  very  good  accuracy  within  3% 
for  a  wide  range  of  problem  configurations.  Hence,  given  the  bounds  //L7pp,  rj^  i  and 
^low'  Tow  i  (3-18)  and  (3.19),  estimators  7upp  and  7iow  of  7*  are  defined  such  that: 


7upp  =  ^(uo,po)  +  I  {(^pp2  “  ??upp2)  +  *(<P,2  “  Vp/^}  ’ 
7low  =  ^(u0 ,  Po)  +  4  {  (^w  -  ,i  -  ^low,i  )  )  > 


(3.25) 


where  7upp  is  derived  by  estimating  the  norms  || ep  ±  £\  \\n  and  ||ep  ±  i£ \  \ \ n  in  the  RHS 
of  (3.17)  by  their  corresponding  upper  bounds  r/^pp,  r/upp  and  where  7iow  is  derived 
similarly  by  using  the  lower  bounds  7^  ,  //|ou,  t. 
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Numerical  experiments  in  Section  3.6  show  that  the  estimators  7avg  and  7|0W  ex¬ 
hibit  poor  effectivity  indices  of  a  factor  10,  for  low  frequencies,  and  a  factor  100,  for 
high  frequencies.  For  low  frequencies,  the  estimator  7upp  gives  the  correct  order  of 
the  magnitude  of  the  error.  However,  for  the  higher  frequency  ranges,  the  effectivity 
index  of  7upp  also  loses  accuracy. 

Toward  deriving  an  alternative  estimator,  the  variational  formulations  in  (3.11) 
governing  the  functions  ei  and  i\  are  recalled,  and  the  sesquilinear  and  hermitian 
properties  of  H(-,  •)  are  employed,  which  yields: 

H(ei  +  £i,v)  =  72(u0,v)  +  72(p0,v),  VvG  V. 

By  comparing  the  explicit  expressions  for  the  primal  and  dual  residual  functional 
in  (3.7),  one  observes  that  H(po,  v)  =  Tl(po,  v).  Hence, 

W(ei  +  ei,v)  =  ft(uo  +  p0,v),  VveF. 

Applying  the  definition  of  the  norm  ||  •  ||^,  as  given  in  (3.13),  leads  to: 

||ei  +  ei||^  =  7e(u0  +  po,ei +£i) 

=  -^(J0(uo  +  Po),ei  +  £i)  +  w2C(j0(uo  +  Po),ei  +£i),  (3.26) 

To  simplify  notations,  in  the  following  treatment,  (uq  +  po)  and  (ei  +  £i)  are  respec¬ 
tively  denoted  as  U  and  e.  From  (3.26)  one  can  establish  the  following  upper  and 
lower  bounds  to  ||£||^: 

||^1(2oU,£)|-|C(joU,£)||  <  }\i\\2n  <  \A(1q\},£)  \  +  |C(j0U,  £)|.  (3.27) 

All  numerical  experiments  in  this  study  have  shown  that  _4(2oU,  e)  <  0  and  C(joU,  e)  > 
0.  This  suggests  that  the  upper  bound  in  the  above  expression  generally  equates  the 
norm  ||e||^.  In  addition,  it  is  observed  that  the  lower  bound  (see  Table  3.2  in  Sec¬ 
tion  3.6.1)  provides  an  accurate  estimator  of  ||e||^  with  effectivity  indices  very  close  to 
1.  Hence,  by  accurately  estimating  |^4(XoCf,e)|  and  \C(jo\J,i)\,  and  by  replacing  these 
terms  in  the  bounds  of  (3.27),  an  accurate  estimate  of  \\i\\n  is  obtained.  Toward  this 
purpose,  the  upper  bound  7+,  as  defined  in  (3.19),  is  recalled: 

m\n  <  4p2  =  ||XoU||3t+llioU||i. 

As  mentioned  previously,  7+  provides  high  accuracies  (within  3%)  for  estimating 
||£||^.  Comparison  of  7+  with  the  upper  bound  in  (3.27),  indicates  that  7+pp  repre¬ 
sents  an  estimate  of  the  upper  bound  in  (3.27)  by  assuming  that  |^4(XoU,e)|  ~  |XnU|:’y 
and  |C(joU,  e)  |  «  ||joU||^.  A  similar  approach  is  proposed  to  estimate  the  lower  bound 
in  (3.27)  by  introducing  £+,  defined  as  follows: 

C+2  =  |l|j0u|l3t-IUoU||?|. 

Numerical  experiments  in  Section  3.6  reveal  that  exhibits  remarkable  accuracy 
within  1%  or  less  for  estimating  ||§i  +  i\\\n  for  a  large  range  of  frequencies  and  prob¬ 
lem  configurations. 
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A  similar  treatment  for  estimating  the  norms  ||§i  —  i]  \\n  and  ||ei  ±  is \\\n  can  be 
applied,  leading  to  the  following  estimates  of  these  norms: 


e±  = 


er 


||J0(u0±po)|l3i-  IUo(uo±Po)|lc 
II  Jo(uo  ±  fpo)  Il3i  -  II  jo(uo  ±  ipo)  lie 


(3.28) 


Subsequently,  these  estimates  are  used  to  propose  an  estimator  7est  of  7*,  defined  as 
follows: 


7est  =  ^(u0,po)  +  I  {(£+2  -  r 2)  +  KCf  -  C2)  }  • 


(3.29) 


This  estimator  is  obtained  by  replacing  the  terms  ||ei  ±  i\ \\n  and  ||§i  ±  ii\  ||^  in  (3.17) 
by  their  respective  estimators,  defined  in  (3.28). 

In  Section  3.6,  numerical  verifications  reveal  that  this  estimator  exhibits  high  ac¬ 
curacy  for  estimating  7*  and  the  modeling  error  itself.  Improving  the  accuracy  of  the 
estimates  of  the  norms  ||  ■  ||^  in  (3.17),  by  using  the  estimates  ^  and  Cf  instead  of  r/hpp 
and  r/dpp ,,  improves  the  overall  accuracy  of  estimating  7*  significantly.  It  is  notewor¬ 
thy  that  for  high  frequency,  the  estimate  7est  still  maintains  good  accuracy. 

Remark  3.5.1  Estimating  |A(XoU,e)|  and  \C(1o\J,e)\  in  (3.27)  zvith  ||J0U||;4  and  ||j0U||^ 
provides  an  estimate  of  the  lozver  bound  in  (3.27),  but  the  estimate  itself  is  not  a  guaranteed 
lozver  bound.  Numerical  experiments  in  Section  3.6  shozv  indeed  that  the  estimates  ^  and 
Cf  are  generally  less  than  ||ei  =t  i\  ||^  and  ||§i  ±  ii\  ||^,  but  in  some  cases  they  can  be  slightly 
larger  (e.g.  see  Tables  3.7  and  3.17). 

Remark  3.5.2  It  is  assumed  that  the  higher  accuracy  of  C±  and  Cf  compared  to  rjfpp  and 
huppi  can  be  explained  by  a  cancellation  of  errors  in  estimates  ||XoU||^  and  |j0U||c  when 
these  are  subtracted  to  estimate  the  lozver  bound  in  (3.27). 


3.5.2  Global  Error  Estimates 

Lemma  3.5.1  Let  (ei,eT)  denote  the  elliptic  representation  defined  in  (3.11)  and  let  (uo,po) 
denote  the  solutions  of  the  coarse  model  (3.4).  Then: 


Clow  —  11^1  Wh  —  Cupp 1 
Clow  —  11^1  llw  X  Cupp  1 
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ivhere: 


Cl 


OW 


Cl 


OW 


P(u0,u0)| 

ll<Mw 

P(Po,Po)l 

llpollw 


Cupp 


Cupp 


■  X0Vu0  +  piv2jo  Uo  •  jo  Uo}  dx, 


{EX0Vpo  -XoVpo  +  puj2jo  po  •  joPo}  dx, 


and  ivhere  the  residuals  functionals  and  deviation  functions  are  given  in  (3.7)  and  (3.8),  re¬ 
spectively. 

Proof:  Only  the  bounds  on  pip  are  proved,  as  the  proof  for  the  bounds  on  pip  is 
similar  by  using  the  approximate  dual  solution  po  instead  of  the  primal  solution  u(l. 
Recalling  the  definition  of  the  norm  ||-  p  and  subsequently  substituting  (3.11)1,  leads 
to: 


|ei||  n  =  sup 
wen 


P(ei,w) 

MP 


=  sup 
wev 


l^-(UQ'W) 

IMP 


The  lower  bound  (jow  follows  quickly  from  this  expression  by  applying  the  definition 
of  the  sup  remum  and  choosing  w  =  uo-  To  prove  the  upper  bound,  recall  the  explicit 
expression  for  the  primal  residual  functional  in  (3.7)  h 


leip  =  sup 
wev 


[  {EX0Vu0  •  Vv  -  puj2jo  u0  •  v}  dx 
Jn 


l|w||  n 


Applying  the  triangle  and  Schwarz  inequalities,  yields: 


eip 


<  sup 
wev 


f  EXoVuo  •  Vwdx 

+ 

/  puj2jo  u0  •  wdx 

Jci 

Jn 

MP 


< 


{EJoVuq  •  X0Vu0  +  pu2j0 u0  •  joUo}  dx. 


Numerical  verifications  of  these  upper  and  lower  bounds  in  Section  3.6  show  that 
they  can  be  used  to  estimate  pip  and  pi  p  within  reasonable  accuracy.  However, 
of  more  interest  are  pop  and  pop.  The  following  lemma  shows  that  these  norms 
are  bounded  from  below  by  their  elliptic  counterparts. 

Lemma  3.5.2  Given  the  solution  pairs  (eoPo)  and  (eipi)  to  (3.9)  and  (3.11),  respectively, 
there  exist  positive  (3 ,  E, uj)  and  Cpfl,  (3 ,  E,  u),  such  that: 
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Proof:  Again,  only  the  first  of  the  two  inequalities  needs  to  be  proved.  Now,  recalling 
the  definition  of  the  norm  ||  •  ||^  and  subsequently  substituting  (3.11) ',  one  obtains: 

„  |^(ei,w)|  |£(e0,w)| 

llei||w  =  sup  — n — n - =  SUP  — n — n - • 

wev  llwllw  wev  IIwIIk 

Substituting  the  definition  of  £>(•,•),  given  in  (2.7),  and  by  applying  the  Schwarz  in¬ 
equality,  leads  to: 


|ei||w<  sup 
wev 


yM(w,w)  +  ca2C(w,w)  +a;H(w,w) 

IHIw 


x  VA(e0,  e0)  +  Lo2C(e0,e0)  +  wP(e0,e0) . 

For  arbitrary  v  G  V,  its  trace  70 v  on  T t  is  in  ,).  Thus: 

^(v,v)  <  ||/3|U‘»(rt)ll7ov|||2(rt)  <  ||/3|U-(rt)||7ov||^i/2(rt). 
The  classical  trace  theorem  for  functions  in  IT1  (FI),  is  used  to  obtain: 
£>(v,v)  <  C(fi)  ||/3||L“(rt)llvllii(o) 

<  (7(0)  ||/3||ioo(rt)  ||E_1||Loo(rt)yl(v,v). 

The  proof  is  completed  by  backsubstituting  this  last  inequality  into  (3.30). 


(3.30) 


In  Section  3.6,  numerical  verifications  are  presented  which  show  that  with  sufficient 
damping,  ||ei||^  and  ||ei||^  are  close  to  ||eo||^  and  pollw-  Hence,  by  estimating  the 
global  norms  of  the  elliptic  representations  of  the  error  functions  within  reasonable 
accuracy,  a  reasonably  reliable  indication  of  the  global  norms  of  the  error  functions 
themselves  is  obtained. 


3.6  Numerical  Experiments 

In  this  section,  several  numerical  verifications  of  the  bounds  and  estimates,  are  pre¬ 
sented.  Sections  3.6.1  through  3.6.3  consider  the  case  of  steady  state  wave  propagation 
for  a  wide  range  of  frequencies.  In  these  sections,  the  effect  of  several  problem  pa¬ 
rameters  on  the  accuracy  of  the  modeling  error  estimators  are  analyzed:  Section  3.6.1 
concentrates  on  the  effect  of  the  impedance  ratio  of  the  elastic  constituents  in  the  mate¬ 
rial,  Section  3.6.2  considers  the  influence  of  damping  boundary  conditions,  and  finally 
Section  3.6.3  shows  the  case  where  the  material  has  two  zones  in  which  the  character¬ 
istic  length  of  the  inhomogeneity  is  different.  A  transient  test  problem  is  given  in 
Section  3.6.4. 


3.6.1  Steady  State  Case  I:  a  Composite  Material  with  Periodic  Microstruc¬ 
ture 

Consider  the  problem  configuration  shown  in  Figure  3.1:  a  beam  with  length  L  is 
clamped  at  its  left  edge  (x  =  0)  and  supported  at  its  right  edge  by  a  damper  with 
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damping  coefficient  /3  =  \J E(L)p(L).  The  beam  has  a  material  microstructure  which 
is  made  out  of  two  elastic  constituents  with  material  properties  {E\,  p{\  and  {E2,  (>2  }, 
that  are  periodically  layered  throughout  the  beam  with  a  constant  layer  thickness  d. 
For  this  test  problem,  the  source  terms  (the  RHS  in  (2.4))  are  characterized  by  an  initial 
displacement  field.  Thus,  Vo(x)  =  0  and  Uo(x)  is  a  symmetric  pulse  located  around 
the  center  point  of  the  beam  xo  with  width  5: 

Uo(x)  =  [x  —  (x0  —  5)] 2  [x  —  (x0  +  (5)] 2 

d  (3.31) 

[1  -  H(x-  (x0  -  5))]  [1  -  H(x  -  (x0  +  (5))] , 

where  H(x  —  a)  denotes  the  Heaviside  function.  The  quantity  of  interest  for  this  nu¬ 
merical  example  is  the  average  stress  on  a  small  domain  S  =  (24, 25)  C  Q: 


To  obtain  the  solution  pairs  (u,  p)  and  (uo,  po)  to  (3.2)  and  (3.4),  respectively,  an  overkill 
computation  is  performed  by  using  approximately  800  quadratic  elements.  The  ap¬ 
proximate  material  model  {iTo,  /Oo}  is  obtained  by  employing  a  standard  classical  asymp¬ 
totic  homogenization  technique  [23].  Hence,  {£o>Po}  are  constant  throughout  the 
beam. 

In  Figures  3.2  through  3.5,  the  solutions  (u,p)  and  (uo,po)  are  shown  for  radial 
frequencies  to  of  500, 1000,  2000,  and  4000  Hz.  In  these  figures,  the  fine  scale  solutions 
(u,p)  are  plotted  as  solid  red  lines  and  the  coarse  scale  solutions  (uo,po)  as  dashed 
green  lines.  These  figures  are  obtained  for  the  case  where  the  beam  consists  out  of 
a  carbon-epoxy  composite;  a  material  that  is  commonly  used  in  engineering  applica¬ 
tions: 


Ei  =  120  GPa,  p\  =  8  g/ cm3,  (Carbon  fiber) 

£'2  =  6  GPa,  P2  =  3  g/cm3,  (Epoxy) 


which  corresponds  to  an  impedance  ratio,  r  = 


lEipi 

E2P2 


,  of  approximately  7.3.  This  value 


is  more  likely  to  trigger  dispersion.  The  corresponding  homogenized  material  prop¬ 
erties  for  the  coarse  scale  problem  (3.4)  are  approximately: 


Eq  =  11.4  GPa,  po  =  5.5  g/cm3. 


A  first  observation  from  Figures  3.2  through  3.5,  is  that  for  the  low  frequencies  of  500 
and  1000  Hz,  the  approximate  solution  uo  ~  u.  In  this  frequency  range,  the  wave 
lengths  are  considerably  larger  than  the  characteristic  length  of  the  inhomogeneity  d. 
Consequently,  the  wave  structure  is  insensitive  to  the  heterogeneous  layers  and  propa¬ 
gates  as  if  the  material  is  homogeneous.  For  these  frequencies,  the  error  in  the  average 
stress  is  entirely  caused  by  the  mismatch  of  E(x)  and  Eo  in  the  domain  of  interest.  For 
higher  frequencies,  however,  the  waves  start  to  notice  the  inhomogeneity.  At  oj  =  2000 
Hz,  the  amplitudes  of  the  two  solutions  already  are  slightly  different.  At  u  =  4000  Hz, 
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L  =  40  m 


(a)  Problem  configuration. 


(b)  Initial  displacement  field 


Figure  3.1:  Steady  state  problem  for  a  composite  material  with  periodic  microstructure. 


19 


0  10  20  30  40 


(a)  Real  part  primal  solution. 


(b)  Imaginary  part  primal  solution. 


(c)  Real  part  dual  solution. 


(d)  Imaginary  part  dual  solution. 


Figure  3.2:  Steady  state  solutions  of  a  periodically  layered  material  with  impedance  ratio  r  =  7.30  and 
for  a  radial  frequency  of  500  Hz,  normalized  by  ||uo||z,“(fi)  and  ||po||L“(f2)- 
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(a)  Real  part  primal  solution.  (b)  Imaginary  part  primal  solution. 


(c)  Real  part  dual  solution. 


(d)  Imaginary  part  dual  solution. 


Figure  3.3:  Steady  state  solutions  of  a  periodically  layered  material  with  impedance  ratio  r  =  7.30  and 
for  a  radial  frequency  of  1000  Hz,  normalized  by  ||uo||L“(n)  an d  l!po||z,“(n)- 
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(a)  Real  part  primal  solution. 


(b)  Imaginary  part  primal  solution. 


Figure  3.4:  Steady  state  solutions  of  a  periodically  layered  material  with  impedance  ratio  r  =  7.30  and 
a  radial  frequency  of  2000  Hz,  normalized  by  ||uo||z,”(q)  and  IIpoHl^o)- 
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(a)  Real  part  primal  solution. 
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(b)  Imaginary  part  primal  solution. 


(d)  Imaginary  part  dual  solution. 


(c)  Real  part  dual  solution. 

Figure  3.5:  Steady  state  solutions  of  a  periodically  layered  material  with  impedance  ratio  r  =  7.30  and 
a  radial  frequency  of  4000  Hz,  normalized  by  ||uo||z,”(q)  and  IIpoHl^o)- 
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the  inhomogeneity  of  the  microstructure  dominates  the  solution.  The  amplitudes  of  u 
and  uo  differ  considerably  and  there  is  a  noticeable  difference  in  phase. 

The  dual  solutions  p  and  po  are  global  functions.  This  is  a  distinctive  difference 
with  the  elastostatic  case  [19,  27,  20,  26],  where  the  dual  solutions  have  very  local  be¬ 
havior,  damping  out  quickly  from  local  responses.  The  ellipticity  of  the  elastostatic 
problem  keeps  the  sensitivity  of  the  solutions  local.  However,  for  the  elastodynamic 
case,  the  hyperbolic  nature  of  the  wave  problem  causes  the  primal  solution  to  be  sen¬ 
sitive  to  global  features.  As  a  consequence,  the  dual  solution  shows  global  behavior. 
This  is  one  of  the  major  complications  in  both  modeling  error  estimation  and  adaptive 
modeling  of  the  wave  problem.  It  requires  a  successful  adaptive  modeling  scheme  to 
be  able  to  perform  nonlocal  adaptation  to  control  the  modeling  error. 

Returning  to  Figures  3.2  through  3.5,  one  sees  that  apart  from  the  amplitude  mis¬ 
matches  between  p  and  po,  the  dual  solutions  show  similar  behavior  as  their  primal 
counterparts.  The  higher  amplitude  mismatch  is  caused  by  the  fact  that  the  force  term 
in  the  dual  problem  has  a  much  larger  amplitude,  due  to  the  presence  of  the  elasticity 
modulus. 

For  a  radial  frequency  c o  =  500  Hz,  Tables  3.1  through  3.3  show  results  on  the 
accuracy  of  the  bounds  and  estimators  derived  in  Sections  3.5.1  and  3.5.2,  To  test  the 
influence  of  the  material  properties,  results  for  three  different  impedance  ratios:  r  = 
1.82,  3.65,  and  7.30  are  presented. 

In  Table  3.1,  the  effectivity  indices  are  listed  for  the  global  bounds  of  Lemma  3.5.1. 
In  the  upper  part,  the  effectivity  indices  with  respect  to  the  global  error  norms  ||ei||ft 
and  pi  ||ft  are  shown.  The  upper  bounds  Cupp  and  CUpP  have  good  effectivity  indices, 
close  to  1.0,  which  improve  as  the  impedance  ratio  increases.  The  lower  bounds  Clow 
and  Ciow  have  poor  accuracy  for  low  impedance  ratios,  but  improve  significantly  as 
the  ratio  increases.  The  lower  part  in  Table  3.1  shows  the  accuracy  with  respect  to  the 
global  error  norms  ||eo lift  and  pollft-  It  is  clear  that  the  global  norms  of  the  elliptic 
representations  (eipi)  are  very  close  to  the  actual  error  norms.  Consequently,  as  the 
effectivity  indices  reveal,  the  upper  bounds  Cupp  and  Cupp  are  close  to  the  error  norms 
as  well  and  represent  accurate  estimators  of  the  global  error  norms,  to  within  2.6  and 
10%  accuracy,  respectively.  In  Table  3.2,  the  effectivity  indices  are  listed  for  the  upper 
and  lower  bounds  riupp(  t)  and  r/^w(.  p  of  Lemma  3.4.1.  Also,  the  effectivity  indices  for 
ajP  are  shown.  These  terms  represent  the  lower  bounds  given  in  (3.27).  The  upper 
bounds  //LCpp,:  t)  are  very  accurate  estimators  within  1.5%.  For  low  r,  the  lower  bounds 
rlk>w(  i)  have  poor  accuracy,  but  they  improve  to  within  5%  accuracy  as  r  increases.  The 
accuracy  of  the  estimators  that  are  proposed  in  Section  3.5.1  and  are  derived  by 
estimating  the  terms  is  remarkable.  Their  accuracy  lies  within  0.5%  and  is  rather 
insensitive  to  the  impedance  ratio  r. 

From  the  results  in  Table  3.2,  one  would  expect  that  the  proposed  error  estimators 
7upp  and  7est  (see  Section  3.5.1)  should  be  accurate  estimators  of  the  modeling  error  in 
the  average  stress.  Table  3.3  lists  the  effectivity  indices  of  these  estimators,  where  the 
effectivity  index  is  defined  by  the  following  ratio: 

,,  .  .  .  ,  | estimator] 

effectivity  index  = - 7 - 

|Q(eo)| 
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r 

1.82 

3.65 

7.30 

Cupp 

1.026 

1.005 

1.026 

Clow 

0.560 

0.709 

0.881 

Cupp 

1.013 

1.005 

1.016 

Clow 

0.618 

0.761 

0.891 

ei  \\h 

0.999 

0.999 

0.999 

|£i \\h 

0.908 

0.902 

0.927 

Cupp 

1.026 

1.005 

1.026 

Clow 

0.560 

0.709 

0.881 

Cupp 

0.921 

0.907 

0.942 

Clow 

0.561 

0.687 

0.826 

Table  3.1:  Effectivity  indices  with  respect  to  |jei||-H  and  ||£i||h  (upper  part),  and  ||eo||-H  and  ||£o||h 
(lower  part),  for  a  periodically  layered  material  and  a  radial  frequency  of  500  Hz. 


r 

1.82 

3.65 

7.30 

^uppGO 

1.013 

1.005 

1.016 

^low(,i) 

0.785 

0.842 

0.914 

ocf 

(.0 

0.998 

0.999 

0.999 

H,i) 

0.995 

1.000 

0.998 

Table  3.2:  Effectivity  indices  of  upper  and  lower  bounds  to  the  terms  coming  from  the  polarization 
formula  expansion  (3.16 ),  for  a  periodically  layered  material,  and  for  radial  frequency  of 
500  Hz. 


r 

1.82 

3.65 

7.30 

|Q(eo)/Q(u)| 

0.807(0) 

0.197(1) 

0.463(1) 

7* 

1.002 

0.977 

1.030 

7uPP 

1.008 

1.145 

3.098 

7low 

4.251 

11.008 

20.050 

7avg 

2.328 

5.108 

8.545 

7est 

1.052 

0.887 

0.976 

Table  3.3:  Relative  error  and  effectivity  indices  of  the  modeling  error  estimators,  for  a  periodically 
layered  material,  and  radial  frequency  of  500  Hz. 
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Indeed,  the  accuracy  for  7est  is  remarkably  good.  The  effectivity  index  for  this  esti¬ 
mator  varies  closely  around  1.0  for  all  impedance  ratios.  The  estimator  7upp  provides 
a  good  estimate  for  low  impedance  ratios  and  becomes  increasingly  inaccurate  as  r 
increases.  As  expected  from  the  results  in  Table  3.2,  the  estimators  7iow  and  7avg  have 
poorer  effectivity  indices  due  to  the  inaccuracy  of  the  lower  bounds  r/Piw(,  .y  Note  that 
7avg  represents  a  decent  estimator  for  low  r  (within  the  order  of  the  error),  but  is  very 
sensitive  to  material  impedance  ratio. 

The  effectivity  indices  for  the  "intermediate"  estimator  7*  ,  which  is  derived  by  elim¬ 
inating  the  §o  terms  from  (3.14)  (see  Section  3.3),  show  that  this  estimator  is  a  good 
intermediate  estimate.  By  ignoring  the  §o  terms,  little  accuracy  is  lost.  The  accuracy 
remains  within  10%. 

Also,  in  Table  3.3,  the  relative  error  is  listed.  Even  at  this  low  frequency,  the  mod¬ 
eling  error  can  be  large:  one  order  higher  than  the  quantity  of  interest  itself.  As  later 
examples  show,  the  error  can  be  a  factor  100  or  1000  higher.  This  is  a  characteristic  fea¬ 
ture  of  the  wave  problem.  The  orders  of  the  modeling  error  can  be  substantially  larger 
than  for  the  elastostatics  case.  This  forms  a  first  indication  that  control  of  the  relative 
modeling  error  to  within  5  or  2%  will  be  computationally  expensive  and,  most  likely, 
not  feasible  in  most  applications.  Tables  3.4  through  3.12  list  the  results  for  higher 
radial  frequencies  of  1000,  2000,  and  4000  Hz.  Comparing  these  results  with  those  we 
discussed  previously  (for  500  Hz),  it  is  observed  from  Tables  3.4  through  3.6  that  the 
effectivity  of  the  global  bounds  Cupp/  Cupp/  Clow/  and  Ci0w  vary  only  slightly  within  an 
accuracy  of  ±5%,  as  the  frequency  increases. 

In  Tables  3.7  through  3.9,  the  effectivity  indices  for  the  bounds  r/^pp(  t  j  and  //^w(  ^ 
show  a  similar  behavior,  with  a  5%  increase  and  decrease,  respectively.  However, 
again  the  effectivity  indices  for  C^{)  behave  remarkably  well.  Only  a  slight  change  of 
accuracy  is  noticed:  in  the  order  of  1%  or  less.  Consequently,  the  accuracy  of  the  corre¬ 
sponding  estimator  7est  remains  good  even  as  the  frequencies  increase  (see  Tables  3.10 
through  3.12).  The  order  of  the  change  of  accuracy  is  much  larger  than  for  the  bounds 


r 

1.82 

3.65 

7.30 

Cupp 

1.024 

1.005 

1.028 

Clow 

0.566 

0.716 

0.880 

Cupp 

1.021 

1.007 

1.022 

Clow 

0.583 

0.758 

0.888 

h  h 

0.998 

0.999 

0.999 

|h \\h 

0.849 

0.892 

0.907 

<%pp 

1.023 

1.004 

1.028 

Clow 

0.565 

0.715 

0.880 

Cupp 

0.868 

0.899 

0.927 

Clow 

0.495 

0.677 

0.805 

Table  3.4:  Effectivity  indices  with  respect  to  |jei||-H  and  ||£i||h  (upper  part),  and  ||g0||^  and  ||£o||h 
(lower  part),  for  a  periodically  layered  material  and  radial  frequency  of  1000  Hz. 
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since  these  terms  estimate  norms  of  large  magnitude.  Hence,  a  small  variation 
in  accuracy  of  is  amplified  by  large  terms,  resulting  in  a  higher  variation  of  the 
accuracy  of  7est.  However,  the  variation  stays  within  an  order  of  1  or  2  with  respect  to 
the  modeling  error,  and  consequently  7est  still  remains  an  accurate  estimator. 


T 

1.82 

3.65 

7.30 

Cupp 

1.025 

1.005 

1.028 

Clow 

0.563 

0.718 

0.880 

Cupp 

1.025 

1.005 

1.027 

Clow 

0.573 

0.724 

0.882 

h  h 

0.997 

0.995 

0.998 

|£i \\h 

0.837 

0.840 

0.887 

Cupp 

1.023 

1.001 

1.027 

Clow 

0.561 

0.715 

0.878 

Cupp 

0.858 

0.845 

0.912 

Clow 

0.480 

0.608 

0.783 

Table  3.5:  Effectivity  indices  with  respect  to  ||§i||w  and  ||ei||«  (upper  part),  and  ||g0||w  and  ||eo||w 
(lower  part),  for  a  periodically  layered  material  and  radial  frequency  of  2000  Hz. 


T 

1.82 

3.65 

7.30 

Cupp 

1.027 

1.005 

1.028 

Clow 

0.562 

0.716 

0.880 

Cupp 

1.027 

1.006 

1.028 

Clow 

0.567 

0.723 

0.881 

h  h 

0.986 

0.995 

0.916 

OH 

0.824 

0.837 

0.883 

Cupp 

1.013 

1.000 

0.942 

Clow 

0.554 

0.713 

0.806 

Cupp 

0.847 

0.843 

0.907 

Clow 

0.468 

0.606 

0.778 

Table  3.6:  Effectivity  indices  with  respect  to  ||@i||«  and  ||£i||-h  (upper  part),  and  ||g0||w  and  ||eo||w 
(lower  part),  for  a  periodically  layered  material  and  radial  frequency  of 4000  Hz. 
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T 

1.82 

3.65 

7.30 

7?± 

'uppGO 

1.021 

1.007 

1.022 

^low  (,i) 

0.635 

0.826 

0.876 

0.987 

1.002 

0.999 

Table  3.7:  Effectivity  indices  of  upper  and  lower  bounds  to  the  terms  coming  from  the  polarization 
formula  expansion  (3.16),  for  a  periodically  layered  material  and  radial  frequency  of 
1000  Hz. 


r 

1.82 

3.65 

7.30 

?7± 

'upp(,0 

1.025 

1.005 

1.027 

^low  (,i) 

0.583 

0.732 

0.843 

P± 

0.984 

0.997 

0.991 

Table  3.8:  Effectivity  indices  of  upper  and  lower  bounds  to  the  terms  coming  from  the  polarization 
formula  expansion  (3.16),  for  a  periodically  layered  material  and  radial  frequency  of 
2000Hz. 


r 

1.82 

3.65 

7.30 

,fupp(,i) 

1.027 

1.006 

1.028 

^low  (,i) 

0.554 

0.727 

0.838 

0.987 

0.998 

0.997 

Table  3.9:  Effectivity  indices  of  upper  and  lower  bounds  to  the  terms  coming  from  the  polarization 
formula  expansion  (3.16),  for  a  periodically  layered  material  and  radial  frequency  of 
4000Hz. 
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r 

1.82 

3.65 

7.30 

|Q(eo)/Q(u)| 

0.860(0) 

0.199(1) 

0.449(1) 

7* 

1.036 

0.998 

0.989 

7uPP 

0.663 

1.014 

2.639 

7low 

16.188 

7.825 

10.623 

7avg 

7.973 

3.975 

4.178 

7est 

1.617 

1.000 

1.138 

Table  3.10:  Relative  error  and  effectivity  indices  of  the  modeling  error  estimators,  for  a  periodically 
layered  material  and  radial  frequency  of  1000  Hz. 


T 

1.82 

3.65 

7.30 

|Q(eo)/Q(u)| 

0.793(0) 

0.219(1) 

0.382(1) 

7* 

1.002 

1.108 

0.732 

7uPP 

1.436 

0.800 

17.993 

7low 

14.257 

46.083 

90.183 

7avg 

6.731 

22.96 

36.096 

7est 

1.063 

1.710 

3.283 

Table  3.11:  Relative  error  and  effectivity  indices  of  the  modeling  error  estimators,  for  a  periodically 
layered  material  and  radial  frequency  of  2000  Hz. 


r 

1.82 

3.65 

7.30 

|Q(eo)/Q(u)| 

0.705(0) 

0.221(1) 

0.425(1) 

7* 

0.858 

1.027 

0.999 

7uPP 

3.981 

1.107 

14.219 

7low 

46.875 

24.932 

78.083 

7avg 

21.509 

12.301 

31.960 

7est 

1.876 

1.114 

0.921 

Table  3.12:  Relative  error  and  effectivity  indices  of  the  modeling  error  estimators,  for  a  periodically 
layered  material  and  radial  frequency  of  4000  Hz. 
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3.6.2  Steady  State  Case  II:  Effect  of  Damping 

In  this  section,  the  effect  of  the  damping  boundary  condition  on  the  accuracy  of  the  es¬ 
timators  is  investigated.  It  is  well  known  that  diminishing  the  damping  in  the  system, 
has  a  destabilizing  effect  on  the  problem  formulation.  To  verify  that  such  a  destabi¬ 
lizing  influence  does  not  affect  the  accuracy  of  the  error  estimators,  a  set  of  numerical 
verifications  is  performed,  using  the  same  problem  configuration  as  in  Section  3.6.1 
(see  also  Figure  3.1),  but  now  the  damping  coefficient  f3  is  decreased  by  multiplying 
this  coefficient  by  a  constant  parameter  (3* . 

In  Table  3.13,  results  for  a  radial  frequency  of  4000  Hz  and  a  set  of  decreasing 
values  of  (3*  are  shown  for  the  effectivity  indices  of  the  global  norm  estimators  of 
Lemma  3.5.1.  The  accuracy  with  respect  to  the  norms  ||ei||^  and  pi||^  appears  to  be 
indifferent  to  the  variation  of  the  damping  coefficient.  However,  the  results  clearly 
indicate  that  as  damping  decreases,  ||ei||^  and  ||ei||w  become  poor  lower  bounds  to 
respectively  ||eo||w  and  polled  especially  ||eo|p.  As  a  direct  result,  the  bounds  Cupp/ 
Ciow/  Cupp/  and  Clow/  become  very  poor  estimators  of  the  global  error  bounds  when 
there  is  very  little  damping  in  the  system. 

Fortunately,  the  accuracy  of  the  bounds  tP, }),  >/^nv( and  is  indifferent  to 
the  variation  of  (3,  as  is  shown  in  Table  3.14.  In  addition.  Table  3.15  illustrates  that  the 
intermediate  estimator  7*,  as  given  in  (3.17),  is  just  slightly  sensitive  to  the  variation 
in  the  damping,  but  still  remains  close  to  1.  Recalling  the  derivations  in  Section  3.5.1 
and  apply  these  two  observations,  one  would  expect  that  the  estimators  7est  and  7upp 
would  maintain  their  accuracy  as  (3*  varies. 

From  Table  3.15,  we  conclude,  however,  that  there  is  a  noticeable  change  in  accu¬ 
racy.  The  accuracy  of  7est  still  remains  acceptable,  yet  has  changed  an  order  of  2.  This 
can  be  explained  by  the  fact  that  the  top  row  in  Table  3.15  indicates  a  dramatic  drop 


(3* 

10° 

lO^1 

10~2 

10~4 

KT8 

Cupp 

1.028 

1.028 

1.028 

1.028 

1.028 

Clow 

0.880 

0.880 

0.880 

0.880 

0.880 

Cupp 

1.028 

1.028 

1.028 

1.028 

1.028 

Clow 

0.881 

0.881 

0.881 

0.881 

0.881 

ei  \\u 

0.916 

0.207 

0.133 

0.132 

0.132 

|h \\h 

0.883 

0.805 

0.719 

0.718 

0.718 

Cupp 

0.942 

0.213 

0.137 

0.135 

0.135 

Clow 

0.806 

0.183 

0.117 

0.116 

0.116 

Cupp 

0.907 

0.827 

0.740 

0.738 

0.738 

Clow 

0.778 

0.709 

0.634 

0.632 

0.632 

Table  3.13:  Effectivity  indices  with  respect  to  ||§i|p  and  ||ei||-H  (upper  part),  and  ||@0||h  and  ||eo||w 
(lower  part),  for  a  periodically  layered  material  and  a  radial  frequency  of 4000  Hz. 
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in  the  relative  modeling  error.  We  now  recall  the  following  equation: 


7*  =  ^(uo,Po)  +  |l|ei  +  ei||w  -  ^||ei  -  ii\\2n 


+^\\e1  +  ii1fn-t-\\el-ie1\\2i. 


(3.17) 


As  shown  in  Section  3.5.1,  the  estimates  jest  and  7upp  are  derived  from  this  expression 
by  estimating  the  norms  in  the  RHS.  These  norms  are  extremely  large  and  as  the  error 
and  7*  become  smaller,  it  will  require  a  higher  accuracy  on  the  estimates  of  the  norms 
to  maintain  the  overall  accuracy  of  the  error  estimators.  Since  the  accuracy  on  the 
norm  estimates  remains  practically  unchanged,  a  loss  in  accuracy  of  the  estimate  of 
the  modeling  error  in  the  average  stress  is  obtained. 

It  is  observed  that  7est  provides  an  excellent  indicator  of  the  modeling  error  for  all 
frequencies. 


(3* 

10° 

KT1 

10~2 

10~4 

10“8 

r ?± 

'upp(,i) 

1.028 

1.028 

1.028 

1.028 

1.028 

^low  (,i) 

0.838 

0.837 

0.837 

0.837 

0.837 

H,i) 

0.997 

0.998 

0.998 

0.998 

0.998 

Table  3.14:  Effectivity  indices  of  upper  and  lower  bounds  to  the  terms  coming  from  the  polarization 
formula  expansion  (3.16),  for  a  periodically  layered  material  and  a  radial  frequency  of 4000 


/ 3 * 

10° 

lO^1 

10~2 

10~4 

10~8 

|Q(eo)| 

|Q(u)l 

0.425(1) 

0.738(0) 

0.414(0) 

0.409(0) 

0.409(0) 

7* 

0.999 

1.197 

1.366 

1.372 

1.372 

7uPP 

14.219 

18.225 

20.826 

20.906 

20.906 

7est 

0.921 

2.072 

2.505 

2.516 

2.516 

Table  3.15:  Relative  error  and  effectivity  indices  of  the  modeling  error  estimators,  for  a  periodically 
layered  material  and  a  radial  frequency  of 4000  Hz. 
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Figure  3.6:  Steady  state  problem  for  a  non-imiformly  layered  material. 


3.6.3  Steady  State  Case  III:  a  Non-Uniformly  Layered  Material 

In  the  previous  examples,  the  material  of  the  beam  has  a  periodic  microstructure.  In 
this  section,  the  more  interesting  nonperiodic  case  is  analyzed.  Consider  the  prob¬ 
lem  configuration  illustrated  in  Figure  3.6:  the  beam  is  constructed  of  a  carbon-epoxy 
composite  material  (see  Section  3.6.1  for  material  properties),  but  now  there  are  two 
different  zones  with  different  layer  thicknesses.  In  the  first  zone,  between  x  =  0  to 
x  =  8,  the  layers  have  thickness  d  =  0.5  m,  whereas  in  the  remaining  part  of  the  beam, 
d  =  0.05  m.  Thus,  between  the  two  zones,  the  layer  thicknesses  differ  a  factor  10  and 
there  should  be  different  dispersion  sensitivities  because  of  this  scale  difference.  In 
the  first  zone,  dispersion  should  be  an  issue  at  lower  frequencies. 

A  damping  boundary  condition  at  the  right  edge  of  the  beam  (/ 3  =  \J E(L)p(L)) 
is  applied  as  before,  but  now  a  nonhomogeneous  Dirichlet  boundary  condition  is  im¬ 
posed  at  the  left  edge,  where  u(0)  =  0.1.  The  source  terms  in  the  RHS  of  (2.4)  are 
all  identically  zero  (no  initial  displacement  or  velocity  field).  For  a  given  frequency 
t o,  the  physical  interpretation  of  this  problem  would  be  that  if  the  edge  at  x  =  0 
is  driven  with  this  frequency  c o,  the  steady  state  solution  of  this  beam  is  given  by 
u (x,t)  =  TZe  {u(x, uj)e'lujt}.  The  quantity  of  interest  for  this  numerical  example  is  the 
average  stress  on  a  small  domain  S  that  is  bordering  with  the  right  edge  of  the  beam: 
S  =  (39,40).  Thus, 


Q(u) 


px=40 
lx= 39 


dx. 


An  overkill  solution  is  computed  by  using  approximately  800  quadratic  elements  to 
compute  the  solution  pairs  (u,p)  and  (uo,po)  to  (3.2)  and  (3.4),  respectively.  Again, 
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the  approximate  material  model  is  determined  by  applying  an  asymptotic  homoge¬ 
nization  technique  and  the  resulting  values  are: 

£'o  =  H-4GPa,  po  =  5.5g/cm3. 

In  Figures  3.7  through  3.10,  the  solutions  (u,p)  and  (uo,po)  are  shown  for  radial  fre¬ 
quencies  uj  of  1000,  2000,  3000,  and  4000  Hz.  In  these  figures,  the  fine  scale  solutions 
(u,p)  are  plotted  as  solid  red  lines  and  the  coarse  scale  solutions  (uo,po)  as  green 
dashed  green  lines.  For  uj  =  1000  and  2000  Hz,  the  wave  lengths  are  significantly  larger 
than  the  dimensions  of  the  inhomogeneity  in  the  material.  Apart  from  a  small  pertur¬ 
bation  in  the  left  zone,  the  exact  solution  u  acts  as  if  it  is  propagating  through  a  ho¬ 
mogeneous  material.  Analogous  to  the  results  for  the  low  frequencies  in  Section  3.6.1, 
the  modeling  error  is  caused  by  the  mismatch  in  E  and  Eq.  For  uj  =  3000  Hz  (see 
Figure  3.9),  the  wave  length  has  reached  such  a  dimension  that  it  is  greatly  disturbed 
by  the  layers  in  the  left  zone  of  the  beam.  Not  only  an  amplitude  mismatch  with  the 
homogenized  solution  is  noticeable,  but  also  a  small  phase  difference.  Whenever  the 
wave  reaches  the  remainder  of  the  beam,  where  the  dimension  of  the  inhomogeneity 
is  much  smaller,  its  wavelength  is  too  large  to  effectively  notice  the  inhomogeneity 
and  propagates  as  if  in  homogeneous  media.  It  is  clear  that  the  dispersion,  created  by 
passing  through  the  left  zone  causes  a  significant  contribution  to  the  error  in  the  aver¬ 
age  stress  in  the  domain  of  interest.  This  example  illustrates  the  global  character  of  the 
wave  problem.  The  accuracy  of  the  solution  in  S  for  this  frequency  is  very  dependent 
on  the  material  features  at  the  other  end  of  the  domain  in  the  left  zone. 

In  Figure  3.10,  the  solutions  for  uj  =  4000  Hz  are  presented  and  here  the  wave 
length  has  reached  a  dimension  such  that  it  is  greatly  affected  by  the  inhomogeneity 
and  attenuates  dramatically  as  it  propagates  through  the  left  zone.  The  frequency  for 
this  wave  is  within  one  of  the  the  so-called  stop-bands  [3,  4]  of  this  inhomogeneity. 
The  solutions  for  the  dual  solution  show  similar  behavior:  for  low  frequencies  there 
is  again  an  amplitude  mismatch  between  p  and  po  (see  Section  3.6.1),  but  they  are 
identical  in  phase.  For  higher  frequencies,  we  again  notice  a  slight  difference  in  phase 
between  p  and  po- 

For  the  set  of  radial  frequencies  c o  of  1000,  2000,  3000,  and  4000  Hz,  Table  3.16 
lists  the  effectivity  indices  for  the  global  bounds  of  Lemma  3.5.1.  The  first  four  rows 
present  the  effectivity  indices  with  respect  to  the  global  error  norms  pi  p  and  pi  p. 
Again,  the  upper  bounds  Cupp  and  Cupp  have  good  effectivity  indices,  close  to  1.  As 
uj  increases,  the  accuracy  deteriorates  slightly  around  2%.  For  Ci0w  and  Clow/  the  ef¬ 
fectivity  indices  are  poorer,  close  to  0.88  and  the  accuracy  diminishes  slightly,  but  still 
within  1%.  The  lower  part  in  Table  3.16  shows  the  accuracy  with  respect  to  the  global 
error  norms  pop  and  pop-  It  is  clear  that  the  global  norms  of  the  elliptic  repre¬ 
sentations  (§i,  ei)  are  very  close  to  the  actual  error  norms  when  the  frequency  is  low, 
but  deviate  as  the  frequency  increases.  Due  to  this  effect,  Cupp  and  Cupp  are  very  ac¬ 
curate  estimators  of  pop  and  pop  for  lower  frequencies,  but  the  corresponding 
accuracies  deteriorate  approximately  15%  as  uj  reaches  4000  Hz.  A  similar  result  is  ob¬ 
served  for  Clow  and  Ciow:  at  uj  =  4000  Hz  their  effectivity  indices  have  reached  a  level 
of  approximately  0.73. 

In  Table  3.17,  the  effectivity  indices  for  the  upper  and  lower  bounds  >l^ppi  t)  and 
Pow(  i)  Lemma  3.4.1  are  presented.  The  upper  bounds  are  very  accurate  estimators 
for  the  low  frequencies,  within  2.0%,  and  lose  2%  accuracy  as  uj  increases.  Again,  the 
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(a)  Real  part  primal  solution. 


(b)  Imaginary  part  primal  solution. 


(c)  Real  part  dual  solution. 


(d)  Imaginary  part  dual  solution. 


Figure  3.7:  Steady  state  solutions  of  a  non-uniformly  layered  material  for  a  radial  frequency  of 
1000 Hz,  normalized  by  ||u0||zJ»(n)  and  ||p0||Loo(n). 


34 


u 

uO 


_ i _ i _ i _ 

0  10  20  30  40 

(a)  Real  part  primal  solution. 
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(b)  Imaginary  part  primal  solution. 
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(c)  Real  part  dual  solution. 


(d)  Imaginary  part  dual  solution. 


Figure  3.8:  Steady  state  solutions  of  a  non-uniformly  layered  material  for  a  radial  frequency  of 
2000Hz ,  normalized  by  ||u0||zJ»(n)  and  ||p0||Loo(n). 
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(a)  Real  part  primal  solution.  (b)  Imaginary  part  primal  solution. 


(c)  Real  part  dual  solution.  (d)  Imaginary  part  dual  solution. 


Figure  3.9:  Steady  state  solutions  of  a  non-uniformly  layered  material  for  a  radial  frequency  of 
3000Hz,  normalized  by  ||u0||zJ»(n)  and  ||p0||Loo(n). 
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(a)  Real  part  primal  solution. 


(b)  Imaginary  part  primal  solution. 


Figure  3.10:  Steady  state  solutions  of  a  non-uniformly  layered  material  for  a  radial  frequency  of 
4000Hz,  normalized  by  ||fi0||z,°°(r2)  and  ||po||L“(fi)- 
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LO  (Hz) 

1000 

2000 

3000 

4000 

Cupp 

1.027 

1.028 

1.029 

1.030 

Clow 

0.880 

0.881 

0.882 

0.882 

Cupp 

1.016 

1.027 

1.029 

1.031 

Clow 

0.891 

0.883 

0.882 

0.884 

ei  \\u 

0.993 

0.990 

0.760 

0.833 

|£l  Wv. 

0.928 

0.882 

0.858 

0.826 

Cupp 

1.021 

1.018 

0.782 

0.859 

Clow 

0.875 

0.872 

0.670 

0.736 

Cupp 

0.943 

0.906 

0.884 

0.852 

Clow 

0.828 

0.779 

0.758 

0.730 

Table  3.16:  Effectivity  indices  with  respect  to  ||@i||«  and  ||ei||-H  (upper  part),  and  ||@0||h  and  ||eo||w 
(lower  part)  for  a  non-uniformly  layered  material. 


to  (Hz) 

1000 

2000 

3000 

4000 

7Cipp(,i) 

1.016 

1.027 

1.029 

1.031 

^low(,i) 

0.913 

0.844 

0.838 

0.840 

1.000 

0.999 

1.000 

1.001 

Table  3.17:  Effectivity  indices  of  upper  and  lower  bounds  to  the  terms  coming  from  the  polarization 
formula  expansion  (3.16),  for  a  non-uniformly  layered  material. 
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lower  bounds  have  poorer  effectivity  indices  than  the  upper  bounds,  and  in  addition 
lose  about  5%  accuracy  as  to  increases.  The  effectivity  indices  for  the  estimators 
once  again  distinguish  themselves  by  both  good  effectivity  indices,  close  to  1,  and  a 
very  small  sensitivity  to  the  frequencies. 

This  directly  affects  the  accuracy  of  the  estimator  7est.  In  Table  3.18,  one  can  see 
that  it  maintains  good  accuracy  for  all  the  frequencies,  with  effectivity  indices  ranging 
from  approximately  1,  for  low  frequencies,  up  to  2.3  for  high  frequencies.  The  slight 
decrease  in  accuracy  is  caused  by  the  accuracy  of  the  intermediate  estimator  7*.  This 
estimator  is  obtained  by  ignoring  C(§o,£i)  (see  (3.14)  and  (3.15)),  which  represents 
an  inertial  feature  within  the  error.  For  higher  frequencies,  this  term  becomes  more 
influential.  However,  the  effectivity  indices  show  that  the  accuracy  is  not  dramatically 
affected  by  eliminating  this  term. 

The  estimators  7avg  and  7iow  again  represent  poor  estimates  and  their  accuracy 
deteriorates  even  further  as  the  frequency  increases.  The  estimator  7upp  shows  a  rea¬ 
sonable  accuracy  for  low  frequencies. 


oj  (Hz) 

1000 

2000 

3000 

4000 

|Q(eo)/Q(u)| 

0.395(1) 

0.471(1) 

0.169(2) 

0.265(5) 

7* 

0.840 

1.021 

2.649 

3.401 

7uPP 

3.709 

11.146 

17.756 

17.945 

7low 

17.104 

60.907 

101.469 

101.419 

7avg 

6.779 

24.899 

41.879 

41.763 

7est 

0.789 

1.176 

2.202 

2.338 

Table  3.18:  Relative  error  and  effectivity  indices  of  the  modeling  error  estimators,  for  a  non-uniformly 
layered  material. 
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3.6.4  Transient  Test  Problem 


In  this  section,  we  present  an  illustrative  example  of  how  the  results  from  the  previ¬ 
ous  steady  state  cases  can  be  applied  to  analyze  transient  problems.  We  consider  the 
carbon-epoxy  composite  beam,  as  shown  in  Figure  3.11.  The  material  microstructure 
is  periodic  uniformly  throughout  the  beam,  where  the  carbon  layers  have  thickness 
d  =  0.05  m,  and  the  epoxy  layers  have  thickness  d  =  0.15  m.  Again,  damping  bound¬ 
ary  conditions  are  prescribed  at  the  right  edge  of  the  beam  and  homogeneous  Dirichlet 
boundary  conditions  at  the  left  edge.  At  t  =  0  s,  the  velocity  field  is  identically  zero,  i.e. 
Vo(x')  =  0,  and  the  initial  displacement  field  Uo(x)  is  prescribed  by  the  pulse  function 
in  (3.31)  with  width  5  =  0.2  m. 

Now,  the  solutions  u {x,t)  and  uo(x,t)  are  computed  for  0  <  t  <  0.05  s,  by  using 
an  overkill  discretization  of  the  frequency  spectrum  of  60  frequencies.  For  each  of 
these  frequencies,  (3.2)  and  (3.4)  are  solved  by  using  an  overkill  mesh  of  800  quadratic 
elements.  The  solutions  in  time  are  subsequently  computed  by  applying  a  discrete 
inverse  Fourier  transformation. 

For  this  specific  example,  the  quantity  of  interest  is  the  average  stress  on  the  small 
subdomain  S  =  (24, 25).  The  approximate  homogenized  material  properties  are  found 
to  be  equal  to: 

Eq  =  7.8  GPa,  po  =  4.2  g/ cm3. 


In  Figure  3.12,  the  solutions  in  time  are  plotted  for  six  consecutive  time  steps.  The  fine 
scale  solution  u(.x,  t)  is  plotted  as  a  solid  red  line  and  the  coarse  scale  solution  uo  0  ,t) 
as  a  green  dashed  line.  Figure  3.12a,  shows  two  major  waves  at  t  =  5  ms  for  both 
solutions:  one  that  propagates  toward  the  left  edge  and  one  toward  the  right  edge.  In 
Figure  3.12b,  the  left-advancing  waves  reach  the  clamped  edge  and  are  about  to  be  re¬ 
flected,  whereas  the  the  other  waves  are  damped  at  the  right  edge.  The  four  following 
graphs  show  the  reflected  waves  propagating  toward  the  right  edge  of  the  beam.  For 
the  later  time  steps  t  =25,  30,  and  35  ms,  the  dispersion  effect  is  very  noticeable.  Com¬ 
pared  to  the  homogenized  solution,  the  initial  peak  of  the  fine  scale  solution  u (x,t) 
is  considerably  smaller.  Also,  the  presence  of  the  trailing  waves  for  u(x,t)  becomes 
apparent.  These  waves  are  much  larger  in  amplitude  than  the  minor  waves  following 
the  initial  pulse  of  the  homogenized  solution  uo(x,t). 

Figure  3.13  shows  the  solutions  in  time  for  the  stresses.  Again,  the  stress  for  the 

d\i(x,  t ) 

fine-scale  problem,  a(x,t)  =  E(x) — y— — ,  is  drawn  as  a  solid  red  line,  whereas  the 

stress  for  the  coarse  scale  problem,  ao(x,t)  =  Eq  — — ,  is  drawn  as  a  dashed  green 
line.  The  solutions  are  shown  for  a  set  of  time  steps  where  the  reflected  waves  are 
passing  through  the  domain  20  <  x  <  30  m.  The  domain  of  interest  S  =  (24, 25)  is  in¬ 
dicated  in  these  graphs  as  well.  A  behavior  similar  to  the  previous  figure  is  observed: 
due  to  the  dispersion  effect,  the  amplitudes  of  the  fine-scale  solution  differ  signifi¬ 
cantly,  and  the  response  is  trailing  behind  with  the  homogenized  solution.  Clearly,  for 
this  example,  a  material  model  based  on  homogenized  material  properties  does  not 
suffice  to  describe  the  wave  problem  accurately. 

Having  computed  Q(§o),  7*,  and  7est;  the  upper  bounds  on  these  quantities  in  the 
time  period  t  €  [0,0.05]  can  be  computed.  This  is  illustrated  on  a  simple  example  of 
an  arbitrary  function  f(x,  t )  for  which  the  Fourier  transform  f(x,  lo)  is  known  at  a  set 
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(a)  Problem  configuration. 


(b)  Initial  displacement  field 


Figure  3.11:  Transient  test  problem  with  initial  displacement  field. 
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(a)  t  =  10  ms. 


0  10  20  30  40 

(b)  t  =  15  ms. 


Figure  3.12:  Transient  solutions  due  to  an  initial  pulse  at  x  =  20  m  with  width  6  =  0.2  m,  normalized 
by  the  maximum  displacement  |uo(x,t)|,  x  €  f2,  t  £  [0, 0.05]. 
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(a)  t  =  30  ms. 


(b)  t  =  31  ms. 


(c)  t  =  32  ms. 


(d)  t  =  33  ms. 


(e)  t  =  34  ms. 


(f)  t  =  35  ms. 


Figure  3.13:  Transient  solutions  for  the  stresses  in  the  domain  of  interest,  normalized  by  the  maximum 

stress  I -Eo -7— (£)£)!/  x  €fl,t  G  [0,0.05]. 
da; 
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4  Adaptive  Modeling  of  the 
Wave  Problem 


The  concept  of  goal-oriented  adaptive  modeling  was  introduced  by  Oden  and  Vema- 
ganti  [19,  20,  26,  27]  for  the  analysis  of  problems  in  elastostatics.  The  basic  philosophy 
of  this  concept  is  to  develop  a  mathematical  model  or  description  out  of  a  hierarchy 
of  available  models,  to  describe  the  occurrence  of  a  physical  event  to  within  a  preset 
solution  accuracy 

In  this  chapter,  this  approach  is  extended  to  the  engineering  problem  of  wave 
propagation  in  heterogeneous  materials.  In  this  case,  the  applications  are  restricted 
to  situations  in  which  the  goal  is  the  average  stress  in  a  small  subdomain  in  the  elastic 
body  In  Section  4.2,  an  adaptive  algorithm  is  proposed  for  the  development  of  the 
material  model  to  control  the  modeling  error.  The  algorithm  employs  error  indicators 
to  decide  where  in  the  material  model  the  adaptation  should  be  implemented.  These 
error  indicators  are  derived  in  Section  4.1.  Finally,  in  Section  4.3,  the  application  of  the 
adaptive  modeling  algorithm  is  demonstrated  for  a  set  of  one-dimensional  examples. 


4.1  Modeling  Error  Indicators 

Recalling  (3.15)  and  applying  the  triangle  and  Schwarz  inequalities,  one  obtains  the 
following: 

|7* I  <  |ft(uo,Po)l  +  |jei||w  pi||w. 

By  applying  Lemma  3.5.1  on  the  norms  in  the  RHS,  gives: 

|7*l<l^(uo,po)|+CupPCUpp-  (4-1) 

Now,  let  the  subdomain  partition  {(-)/,  }  of  f!  be  defined  as  follows: 

Q  =  int  ^  [J  0  ,  Qi  IT  Qj  =  0,  Vi  /  j.  (4.2) 

Then  (4.1)  suggests  that  the  restrictions  to  { 0 /, }  of  each  of  the  terms  in  the  RHS  give 
an  indication  of  the  contribution  of  the  subdomains  to  17*1.  Since  7*  forms  a  reliable 
estimate  of  the  modeling  error  (see  Section  3.6),  these  indicators  can  be  extrapolated 
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to  estimating  the  modeling  error  itself.  Thus,  the  following  quantities  {£&},  defined 
on  {©fc},  are  proposed  as  error  indicators  of  the  modeling  error: 


8k  =  |7^(uo,po)|  +  Ck  Ck, 


where 


ftfc(u0,Po)  =  /  {EXoVuo  •  Vp0  —  puj2 jo  uo  •  po}  dx, 


(4.3) 


C fc  =  J  J_  {ETqVuo  •  ToVuq  +  p<jj2jo  Uq  •  jo  uo}  dx, 


Ck  =  J {EToVpo  •  J0Vp0  +  pw2j0po  •  JoPojdx, 


where  the  deviation  tensor  Xq  and  function  jo  are  given  in  (3.8). 


4.2  The  Adaptive  Algorithm 

In  this  section,  a  detailed  description  of  the  adaptive  modeling  algorithm,  proposed 
for  the  analysis  of  wave  propagation  through  heterogeneous  materials,  is  given.  The 
occurrence  of  geometric  dispersion  [3,  4]  is  a  major  complication  in  modeling  elastic 
wave  propagation.  Contrary  to  the  elastostatic  case,  the  sensitivity  of  the  solutions  is 
very  global,  as  the  numerical  examples  in  Section  3.6  show. 

Since  dispersion  is  determined  by  the  ratio  between  the  wave-length  and  the  di¬ 
mension  of  the  inhomogeneity,  it  follows  directly  that  there  is  a  great  dependence  on 
the  wave  frequency.  Therefore,  the  choice  for  a  frequency  domain  analysis  seems  to  be 
reasonable.  The  basic  idea  is  that  by  decomposing  the  wave  into  a  spectrum  of  steady 
state  waves  with  frequency  u,  those  frequencies  are  identified  that  are  sensitive  to  the 
material  inhomogeneity.  Subsequently,  the  adaptation  of  the  material  model  can  be 
applied  for  these  frequencies. 

In  the  following,  a  step  by  step  account  of  the  adaptive  modeling  algorithm  is 
given,  which  is  illustrated  schematically  in  Figure  4.1: 

Step  1.  The  domain  O  is  partitioned  into  subdomains  {©fc},  as  given  by  (4.2),  where 
one  of  the  subdomains  is  S,  the  area  on  which  the  quantity  of  interest  is  defined 
(see  (3.1)).  The  set  {©fc}  is  chosen  such  that  each  subdomain  contains  inhomo¬ 
geneities  of  similar  dimensions. 

In  this  study,  the  coarse  model  {Xo,  po}/  is  chosen  as  the  homogenized  (constant) 
material  properties  that  are  obtained  by  applying  a  classical  homogenization 
technique  [23].  In  general  applications,  the  user  has  the  freedom  to  choose  any 
abstract  approximate  material  to  start  the  algorithm. 

Step  2.  The  wave  problem  is  decomposed  into  a  set  of  steady  state  waves,  with  fre¬ 
quencies  {ujn  },  where  n  =  1, 2, . . .  TV.  For  the  analysis  of  one  steady  state  wave, 
of  course  N  =  1,  whereas  for  a  wave  with  a  discrete  frequency  spectrum  we  have 
N  >  1.  For  the  analysis  of  transient  waves  with  continuous  frequency  spectra. 
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Figure  4.1:  Flozv  diagram  of  the  adaptive  modeling  algorithm. 
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an  overkill  discretization  of  the  frequency  spectrum  is  used.  The  overkill  dis¬ 
cretization  is  assumed  to  deliver  the  desired  accuracy 

Step  3.  For  every  frequency  uin,  the  solution  pair  (uo,po)  to  (3.4)  is  computed.  Since 
the  numerical  experiments  in  Section  3.6  show  that  the  estimator  7est  (see  Sec¬ 
tion  3.5.1)  has  highest  accuracy,  this  estimator  is  chosen  to  assess  the  modeling 
error  and,  consequently,  is  computed  for  every  frequency  as  well. 

Step  4.  By  using  the  {7est}  for  all  the  frequencies,  the  upper  bound  in  time  on  Q(eo)|, 
the  error  in  the  average  stress  on  S,  is  computed  by  applying  (3.33).  In  the  case 
where  N  =  1,  only  1 7est |  needs  to  be  computed. 

Step  5.  The  estimate  of  the  error  bound  is  compared  with  a  user-set  tolerance.  If  the 
estimate  exceeds  the  error  tolerance,  the  analysis  proceeds  to  Step  6,  if  not,  the 
analysis  ends. 

Step  6.  A  search  through  the  frequency  spectrum  is  performed  to  identify  the  frequency 
wcrit  that  has  the  highest  contribution  to  the  error  bound  of  Step  4. 

Step  7.  Subsequently,  the  error  indicators  {£/,.},  as  given  in  (4.3),  are  computed  for  all 
the  subdomains  {©&}. 


Step  8.  Next,  for  tucrit,  the  material  model  is  adapted  in  the  subdomain  0/,.  that  has  the 

highest  relative  error  indicator  — — .  In  this  study,  the  adaptation  implies  the 

|0fc| 

inclusion  of  the  exact  material  properties  {E,p}  in  0/,.  Of  course,  in  general 
applications,  one  can  choose  different  adaptation  techniques,  e.g.  by  using  a 
hierarchy  of  material  properties.  It  is  emphasized  that  at  the  end  of  the  process, 
the  approximate  material  model  is  generally  different  for  different  frequencies. 


Step  9.  Finally,  for  u,’cnt  only,  the  solution  pair  (uq,  po)  to  (3.4)  is  recomputed  by  employ¬ 
ing  the  updated  material  model  {Eq(x): po(x)},  and  7est  is  recomputed  for  this 
frequency.  Then  Step  5  is  revisited. 


4.3  Numerical  Examples 

In  this  section,  a  set  of  numerical  examples  is  presented  to  demonstrate  the  application 
of  the  adaptive  algorithm  proposed  in  Section  4.2.  The  application  to  steady  states 
wave  of  one  frequency  is  given  in  Section  4.3.1,  whereas  Section  4.3.2  treats  transient 
waves  that  are  bandlimited. 


4.3.1  Steady  State  Waves 

The  problem  configuration  considered  is  shown  in  Figure  4.2:  a  beam  with  3  zones 
with  different  dimensions  of  inhomogeneity.  As  in  Section  3.6,  the  material  is  made 
out  of  two  constituents:  carbon  and  epoxy  (see  Section  3.6  for  material  properties).  The 
beam  has  a  damping  boundary  condition  at  its  right  edge,  and  a  nonhomogeneous 
Dirichlet  boundary  condition  u(0)  =  0.1  at  the  left  edge. 
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The  beam  is  partitioned  into  7  subdomains,  such  that: 


0! 

=  (0,8), 

(contains  zone  I) 

02 

=  (8,16), 

(in  zone  II) 

©3 

=  (16,24), 

(in  zone  II) 

04 

=  (24,32), 

(in  zone  II) 

05 

=  (32,36), 

(in  zone  III) 

©6 

=  (36.0,36.1), 

(domain  of  interest) 

07 

=  (36.1,40). 

(in  zone  III) 

The  domain  of  interest  S,  on  which  the  average  stress  is  computed,  is  represented  by 
06-  This  small  subdomain  consists  of  only  two  layers  and  is  a  representative  cell  for 
the  inhomogeneity  in  zone  III.  The  approximate  material  properties  at  the  beginning 
of  the  adaptive  process,  {Eq.  p0),  are  obtained  by  applying  the  classical  asymptotic 
homogenization  technique  and  using  the  domain  of  interest  S  as  the  representative 
unit  cell  of  the  entire  beam.  Thus, 

l?o  =  11.4GPa,  po  =  5.5  g/cm3. 

The  adaptive  algorithm  is  first  applied  to  the  case  in  which  the  radial  frequency  is  low, 
i.e.  lu  =  200  Hz.  In  Table  4.1  a  data  log  is  shown  for  every  iteration  step  in  the  algo¬ 
rithm.  It  lists  the  relative  error,  the  effectivity  index  for  the  estimator  7est,  and  the  list 
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of  the  distribution  of  the  relative  error  indicators.  The  subdomain  that  has  the  highest 
relative  error  indicator  is  highlighted  red  and  for  the  next  iteration  step,  the  material 
properties  are  adapted  in  these  domains  by  taking  the  exact  material  properties  (see 
Section  4.2).  From  the  results  for  the  relative  errors,  it  is  observed  that  the  problem 
behavior  is  still  rather  local  for  this  low  frequency.  A  drop  from  a  relative  error  of  4.5 
to  0.03  is  observed.  This  significant  increase  in  accuracy  is  achieved  by  just  including 
the  exact  material  properties  in  the  domain  of  interest.  For  iteration  1  and  3,  it  is  ob¬ 
served  that  the  effectivity  index  of  7est  exhibits  good  accuracy.  However,  for  iteration 

2,  a  poorer  effectivity  of  approximately  3  is  observed.  This  behavior  is  observed  when 
the  relative  error  drops  below  10%.  The  effectivity  indices  range  from  1  to  3  (though 
an  effectivity  index  of  3  is  rarely  seen).  This  is  due  to  the  fact  that,  though  the  error 
is  small,  the  estimator  involves  the  difference  between  terms  ||ei  ±  (i)i 1 \\n  (see  (3.29)) 
that  still  remain  very  large.  Since  the  accuracy  of  the  estimates  of  these  terms, 
does  not  improve  proportionally  (see  also  results  in  Tables  3.14  and  3.15),  accuracy  is 
lost  in  the  overall  estimator  7est.  Still,  the  estimator  is  of  the  right  order.  After  iteration 

3,  the  adaptive  process  is  stopped,  as  a  very  high  accuracy,  of  less  than  1%,  is  achieved. 

For  the  three  iterations.  Figure  4.3  shows  the  real  and  imaginary  parts  of  the  solu¬ 
tions,  where  the  fine  scale  solution  u  is  plotted  as  a  solid  red  line  and  the  approximate 
solution  uo  as  a  green  dashed  green  line.  The  graphs  on  the  left  show  that  the  coarse 
material  model  for  Eq  (the  profiles  for  fx,  are  similar)  changes  with  each  iteration  step. 


Iteration 

|Q(e0)|  |7est| 

|Q(u)l  |Q(eo)| 

£ 

Error  Indicators  — — - 

|Bfc| 

1 

0.452(1)  0.98 

1  -  0.910(9) 

2  -  0.303(9) 

3  -  0.704(9) 

4  -  0.810(9) 

5  -  0.208(9) 
6-0.128(12) 
7-0.200(9) 

2 

0.298(— 1)  3.08 

1  -  0.165(9) 

2  -  0.548(8) 

3- 0.127(9) 

4- 0.147(9) 

5  -  0.357(8) 
7-0.368(8) 

3 

0.335(— 2)  1.31 

2  -  0.558(8) 

3  -  0.128(9) 
4-0.154(9) 

5  -  0.386(8) 

7  -  0.329(8) 

Table  4.1:  Data  review  of  the  adaptive  modeling  process  for  the  steady  state  wave  ivitli  a  radial 
frequency  of  200  Hz. 
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Iteration  1: 


Ep(x) 

l|£|U~(n) 


-  Real  parts 


-  Imaginary  parts 


Iteration  2: 


E0{x) 

ll-BIU‘»(n) 


-  Real  parts 


-  Imaginary  parts 


Iteration  3: 


E0(x) 


-  Real  parts 


||-E|U°°(n) 


10  20  30  40 

-  Imaginary  parts 


Figure  4.3:  Snapshot  solutions  of  the  adaptive  modeling  analysis  for  the  steady  state  wave,  ivith 
u>=200  Hz,  normalized  by  ||uo||_L«=(n)- 
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In  Table  4.2,  the  adaptation  results  are  listed  for  the  more  interesting  case  in  which 
the  frequency  is  high,  u  =  3000  Hz.  The  relative  error  again  decreases  with  each  it¬ 
eration  and  between  some  steps  a  factor  10  difference  is  observed.  For  iteration  1 
through  6,  the  effectivity  indices  for  7est  exhibit  good  accuracy,  with  values  close  to  1. 
For  the  last  iteration  step,  in  which  the  relative  error  is  <10%,  a  loss  in  accuracy  for 
7est  is  observed,  but  the  estimator  is  still  of  the  right  order.  The  explanation  for  the 
loss  in  accuracy  is  identical  to  the  case  in  which  a ;  =  200  Hz,  iteration  2  (see  previous 
remarks  on  Table  4.1). 

Figure  4.4  displays  the  coarse  material  model,  and  the  real  and  imaginary  parts 
of  the  solutions  u  and  Uo  for  iteration  steps  1,  4,  and  7.  For  iteration  4  and  7,  Fig¬ 
ure  4.4  shows  that  the  adaptive  process  has  incorporated  the  exact  material  properties 
in  some  subdomains  into  the  coarse  model. 


Iteration 

|Q(eo)| 

I  'Test  | 

|2(u)| 

|Q(eo)| 

Bfc 

1 

0.152(2) 

1.44 

1-0.521(12) 

2  -  0.468(12) 
3-0.481(12) 

4  -  0.503(12) 

5  -  0.445(12) 

6  -  0.424(13) 

7  -  0.577(12) 

2 

0.354(1) 

0.85 

1-0.928(11) 

2  -  0.835(11) 

3- 0.858(11) 

4- 0.896(11) 

5  -  0.794(11) 
7-0.106(12) 

3 

0.353(1) 

1.41 

1-0.885(11) 

2  -  0.796(11) 

3- 0.820(11) 

4- 0.854(11) 

5- 0.758(11) 

4 

0.324(0) 

1.19 

2  -  0.244(11) 

3  -  0.246(11) 

4- 0.260(11) 

5- 0.231(11) 

5 

0.278(0) 

0.94 

2- 0.196(11) 

3- 0.197(11) 

5  -  0.222(11) 

6 

0.253(0) 

0.81 

2-0.228(11) 

3  -  0.229(11) 

7 

0.889(— 1) 

2.33 

Table  4.2:  Data  review  of  the  adaptive  modeling  process  for  the  steady  state  wave  ivitli  a  radial 
frequency  of  3000  Hz. 
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Iteration  1: 


Ep(x) 

ll-E|U“(n) 


-  Real  parts 


-  Imaginary  parts 


E0(x) 


10 


Iteration  4: 


20  30 

Ep(x) 

'  l|s|U-(n)- 


40 


-  Real  parts 


-  Imaginary  parts 


Iteration  7:  - 


Ep(x) 


-  Real  parts 


-  Imaginary  parts 


Figure  4.4:  Snapshot  solutions  of  the  adaptive  modeling  analysis  for  the  steady  state  wave,  with 
u>=3000  Hz. 
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Note  that  the  scale  used  in  these  graphical  representations  is  such  that  the  highly  os¬ 
cillatory  material  properties  of  adapted  subdomains  in  zone  II  and  III  appear  to  be 
constant.  For  this  frequency,  capturing  dispersion  plays  an  important  role  in  achiev¬ 
ing  an  accurate  solution.  At  the  first  iteration,  with  all  homogeneous  approximate 
material  properties,  there  is  a  significant  mismatch  between  the  exact  and  approxi¬ 
mate  solution.  The  results  for  iteration  4,  reveal  that  the  sensitivity  of  the  wave  to  the 
heterogeneity  in  zone  I  (or  ©i)  is  the  main  source  of  error  at  the  first  iteration. 

The  step  between  iteration  3  and  4  clearly  shows  the  nonlocal  behavior  of  the  adap¬ 
tation  process.  After  the  material  model  in  ©i,  far  away  from  the  domain  of  interest 
S,  is  adapted,  a  considerable  increase  of  accuracy  of  the  approximate  solution  is  ob¬ 
served. 

Finally,  for  a  radial  frequency  of  4000  Hz,  Table  4.3  and  Figure  4.5  display  simi¬ 
lar  results  as  seen  for  uj  =  3000  Hz.  Again,  7est  exhibits  good  accuracy  and  nonlocal 
model  adaptation  appears  between  iteration  3  and  4,  and  iteration  5  and  6.  An  ex¬ 
tremely  large  range  of  the  relative  error  is  noted.  For  iteration  1,  a  relative  error  of 
approximately  3,000,000%  is  observed,  which  decreases  continuously  as  the  adaptive 
process  is  performed  until  it  reaches  a  value  of  approximately  10%  at  iteration  7.  For 
analyses  in  which  the  computational  cost  is  a  dominant  criterion  over  solution  accu¬ 
racy,  this  final  value  of  10%  suggests  that  error  tolerances  in  the  range  of  5%  may  not 
be  practical. 
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Iteration 

|Q(e0)|  |7est| 

|2(u)|  |Q(e0)| 

£k 

Error  Indicators  — — r 

\wk\ 

1 

0.292(5)  1.51 

1  -  0.550(12) 

2  -  0.539(12) 

3  -  0.538(12) 

4  -  0.537(12) 

5  -  0.505(12) 

6  -  0.506(13) 

7  -  0.567(12) 

2 

0.523(4)  0.84 

1-0.971(11) 

2  -  0.953(11) 

3- 0.951(11) 

4- 0.950(11) 

5- 0.839(11) 
7-0.110(12) 

3 

0.479(1)  1.39 

1-0.800(11) 

2  -  0.784(11) 

3- 0.781(11) 

4- 0.781(11) 

5- 0.739(11) 

4 

0.609(0)  1.03 

2  -  0.254(8) 

3  -  0.253(8) 

4  -  0.254(8) 

5  -  0.263(8) 

5 

0.304(0)  0.65 

2-0.196(8) 

3  -  0.195(8) 
4-0.195(8) 

6 

0.278(0)  0.81 

3-0.172(8) 

4  -  0.173(8) 

7 

0.127(0)  0.99 

Table  4.3:  Data  review  of  the  adaptive  modeling  process  for  the  steady  state  wave  with  a  radial 
frequency  of 4000  Hz. 
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E0(x)  - 
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u  - 

ful  1  i  ■'  i 
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11  1  *  11  n  *  *i  h  jfi 
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11  y 

il  V  if  V  l!  <  V  v  Si 

;  j!  ;!  !  v  V'  SI  •;  If  l 

'<  l1  i 

ii  ii  ii  ii !«'  ii  ij  ii  ii  ii  si  s 

IJ  j;  i!  I,'  ;l  S|  ;!  ;l  ;|'  SI  i;  ;i 

J  <1  SI  il  Si  J  V  If  V  SI  v  i 

0  10  20  30  40  0  10  20  30  40  0  10  20  30  40 


Iteration  1:  - 


Eq(x) 

ll-E|U“(n) 


-  Real  parts 


-  Imaginary  parts 


o 


10 


E0(x) 


20  30  40 


0  10  20  30  40 


Iteration  3:  - 


Ep(x) 

||-E|[r,~(n) 


-  Real  parts 


-  Imaginary  parts 


Iteration  4:  - 


Ep(x) 


-  Real  parts 


-  Imaginary  parts 


Figure  4.5:  Snapshot  solutions  of  the  adaptive  modeling  analysis  for  the  steady  state  wave,  with 
u=4000  Hz. 
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4.3.2  Transient  Bandlimited  Waves 


The  problem  configuration  is  shown  in  Figure  4.6:  a  beam  with  4  zones  with  different 
dimensions  of  inhomogeneity.  Again,  the  material  is  made  out  of  layers  of  carbon  and 
epoxy  (see  Section  3.6  for  material  properties).  The  beam  has  a  damping  boundary 
condition  at  its  right  edge,  and  a  driven  displacement  (Dirichlet)  boundary  condition 
at  the  left  edge,  whose  frequency  spectrum  is  bandlimited  with  a  maximum  frequency 

^max/ 

for  |iu|  <  Wmax, 
for  |(u|  >  LUmax- 

The  domain  is  partitioned  into  9  subdomains,  such  that: 


0! 

=  (0,4), 

(contains  zone  I) 

02 

=  (4,10), 

(in  zone  II) 

©3 

=  (10,16), 

(in  zone  II) 

04 

=  (16,22), 

(in  zone  II) 

©5 

=  (22,25), 

(contains  zone  III) 

©6 

=  (25,30), 

(in  zone  IV) 

07 

=  (30,31), 

(domain  of  interest) 

©8 

=  (31,35), 

(in  zone  IV) 

©9 

=  (35,40). 

(in  zone  IV) 

The  domain  of  interest  S,  on  which  the  average  stress  is  sought,  is  represented  by 
@7.  The  approximate  material  properties  at  the  beginning  of  the  adaptive  process  are 
again  obtained  by  applying  the  classical  asymptotic  homogenization  technique  and 
using  the  domain  of  interest  S  as  the  representative  unit  cell  of  the  entire  beam.  Thus, 

Eq  =  14.2  GPa,  po  =  6-Og/cm3. 

The  case  in  which  ccmax  is  relatively  low,  ccmax  =  1000  Hz,  is  considered.  The  adap¬ 
tive  algorithm  is  applied  by  using  a  discretization  of  the  frequency  spectrum  into  10 
frequencies.  The  radial  frequency  un  of  a  wave  with  frequency  number  n,  is  given  by: 

2irn 

u)n  =  —,  n  =  0,1,...,  9 

where  the  end  time  of  the  analysis  T  is  57  ms.  In  Table  4.4,  an  overview  of  the  results 
after  the  first  iteration  step  is  shown,  for  9  of  the  10  frequencies.  Here,  the  approximate 
material  model  consists  of  homogenized  material  properties.  The  results  for  n  =  0  are 
not  listed  as  this  mode  does  not  have  a  contribution  to  the  quantity  of  interest  (a  zero 
estimate  is  therefore  100%  accurate). 
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By  applying  (3.33)  upper  bounds  in  time  are  obtained  on  the  error  in  the  quan¬ 
tity  of  interest  and  the  quantity  of  interest  itself.  For  this  specific  case,  the  ratio  of 
the  upper  bound  on  the  error  and  the  bound  on  the  quantity  of  interest  is  9.1  after 
the  first  iteration  step,  which  suggests  high  relative  error  for  the  time  frame  [0  ,T\. 


frequency  no. 

c o  (Hz) 

|Q(eo)| 

l7est| 

iswi 

IS(*o)l 

1 

111 

4.94 

1.07 

2 

222 

7.56 

1.35 

3 

333 

6.00 

1.55 

4 

444 

8.70 

1.96 

5 

556 

6.73 

1.25 

6 

667 

13.3 

1.40 

7 

778 

11.8 

1.06 

8 

889 

10.9 

1.67 

9 

1000 

4.34 

1.30 

Table  4.4:  Relative  errors  and  effectivity  indices  for  the  estimator  yest,  after  the  first  iteration  step, 
when  Umax  =  1000  Hz. 
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By  using  the  estimates  {7est}/  an  upper  bound  in  time  on  the  error  estimate  is  obtained. 
After  the  first  iteration  step  of  the  of  the  adaptive  loop,  the  ratio  of  the  bounds  on 
the  error  estimate  and  the  error  itself  is  1.38,  which  represents  good  accuracy  of  the 
estimator.  The  adaptive  modeling  algorithm  is  applied  until  the  ratio  of  the  upper 
bounds  on  the  error  and  quantity  of  interest  reaches  a  value  below  0.4.  In  Table  4.5, 
an  overview  of  the  entire  adaptive  iteration  procedure  is  recorded.  For  every  iteration 
step,  the  critical  frequency  number,  the  subdomain  with  the  highest  contribution  to 
the  error,  the  ratio  of  the  upper  bounds  on  the  error  and  the  quantity  of  interest,  and 
the  effectivity  index  of  the  error  estimator  are  listed.  Note  that  the  latter  is  defined  as 
follows: 

„  .  .  T  .  upper  bound  on  the  error  estimator 

bff ectivity  Index  =  — LJ - - - - - - - . 

upper  bound  on  the  error 

Initially,  it  is  the  domain  of  interest  @7  where  the  adaptation  is  applied  for  the  first  7  it¬ 
erations.  The  adaptation  process  subsequently  concentrates  on  other  subdomains.  It  is 
observed  that  especially  for  the  higher  frequencies  the  material  model  is  adapted.  Ap¬ 
parently,  these  waves  are  most  sensitive  to  the  inhomogeneity  in  the  material.  Overall, 
the  estimate  of  the  upper  bound  on  the  error  shows  good  accuracy.  Toward  the  end  of 
the  adaptive  process,  the  effectivity  index  is  pushed  toward  2.  This  is  caused  by  the 
effect  observed  previously:  in  those  cases,  there  are  waves  for  which  the  relative  error 
is  below  10%,  and  for  which  the  effectivity  index  of  the  estimator  7est  moves  toward 
an  effectivity  of  2  (see  remarks  on  Tables  4.1  and  4.2  in  Section  4.3.1).  As  a  result,  the 
overall  effectivity  index  shifts  to  slightly  higher  values. 

The  solutions  in  time  are  computed  by  applying  the  discrete  inverse  Fourier  trans¬ 
formation,  as  given  in  (3.32).  The  results  for  the  computed  stresses  at  times  t  =14, 16, 
and  18  ms  are  plotted  in  Figure  4.7.  In  these  graphs  the  exact  stress  is  drawn  as  a 
solid  red  line.  In  the  graphs  in  the  left  column,  the  approximate  solution,  as  obtained 
after  the  first  iteration  (by  using  the  homogenized  material  properties),  is  drawn  as 
a  dashed  blue  line.  The  approximate  solution,  obtained  at  the  end  of  the  adaptive 
modeling  process,  is  drawn  as  a  dashed  green  line  in  the  graphs  in  the  right  column. 
At  t  =  14  and  18  ms,  it  is  observed  that  the  mismatch  in  the  homogenized  and  ex¬ 
act  solution  is  large  on  the  domain  of  interest  S  (also  indicated  in  these  graphs).  The 
"adapted"  solution,  however,  is  very  accurate  and  follows  the  exact  solution  quite 
nicely. 

Results  where  the  maximum  frequency  is  higher,  namely:  umax  =  2000,  3000,  and 
4000  Hz,  are  also  obtained.  For  these  cases,  the  results  after  the  first  and  the  last  itera¬ 
tion  of  the  adaptive  process  are  listed  in  respectively  Tables  4.6  and  4.7.  At  the  end  of 
the  first  iteration,  the  ratios  of  the  upper  bounds  on  the  error  and  the  quantity  of  in¬ 
terest  are  large,  close  to  7  for  all  three  cases,  which  suggest  large  relative  errors  for  the 
time  frame  [0,  T].  The  estimate  of  the  upper  bound  on  the  error  appears  to  be  accurate 
and  the  effectivity  index  varies  between  1.1  to  1.3.  At  the  end  of  the  adaptive  process, 
the  relative  error  has  been  reduced  significantly,  since  the  ratio  of  the  upper  bounds 
on  the  error  and  the  quantity  of  interest  is  low,  close  to  0.3,  for  all  three  frequencies. 
The  effectivity  index  of  the  estimator  has  shifted  to  slightly  higher  values  closer  to  2. 
For  the  cases  with  higher  1 omax,  more  iterations  are  needed  to  obtain  the  requited  ac¬ 
curacy.  This  is  caused  by  the  fact  that  the  inhomogeneity  in  this  problem  has  a  small 
scale.  Only  waves  with  small  wave  length  (thus,  high  frequencies)  will  be  sensitive  to 
the  inhomogeneity. 
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Iteration 

No. 

Critical  Freq. 
No. 

@crit 

|C(eo)l 

|Q(u)| 

Eff .  Index  of  estimate 
on  upper  bound 

1 

7 

7 

7.09 

1.45 

2 

6 

7 

5.41 

1.43 

3 

8 

7 

3.08 

1.25 

4 

9 

7 

2.90 

1.41 

5 

5 

7 

2.33 

1.25 

6 

3 

7 

1.81 

1.17 

7 

4 

7 

1.40 

0.94 

8 

6 

2 

1.32 

1.04 

9 

4 

7 

1.05 

1.56 

10 

7 

9 

1.01 

1.04 

11 

7 

4 

0.94 

1.05 

12 

8 

4 

0.86 

1.00 

13 

6 

9 

0.81 

1.03 

14 

7 

3 

0.75 

1.11 

15 

6 

6 

0.72 

1.12 

16 

8 

3 

0.66 

1.18 

17 

6 

3 

0.64 

1.18 

18 

6 

8 

0.63 

1.12 

19 

7 

2 

0.58 

1.18 

20 

1 

7 

0.48 

1.21 

21 

6 

5 

0.45 

1.26 

22 

5 

6 

0.44 

1.50 

23 

5 

4 

0.44 

1.57 

24 

5 

2 

0.43 

1.69 

25 

8 

6 

0.39 

1.79 

Table  4.5:  Summary  on  iterative  adaption  process  for  unuix=1000  Hz. 
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^max  (Hz) 

Number  of 
frequencies 

|G(eo)| 

|Q(u)| 

Eff.  Index  for  estimate 
on  the  upper  bound 

2000 

19 

7.35 

1.31 

3000 

28 

6.99 

1.10 

4000 

37 

7.46 

1.18 

Table  4.6:  Ratios  of  the  upper  bounds  on  the  error  and  the  quantity  of  interest,  and  effectivity  indices 
for  the  upper  bounds  on  the  error  estimator,  after  the  first  iteration  step. 


^max  (Hz) 

Number  of 
iterations 

|Q(eo)| 

|G(u)| 

Eff.  Index  for  estimate 
on  the  upper  bound 

2000 

47 

0.35 

1.32 

3000 

93 

0.33 

1.64 

4000 

128 

0.32 

2.33 

Table  4.7:  Ratios  of  the  upper  bounds  on  the  error  and  the  quantity  of  interest,  and  effectivity  indices 
for  the  upper  bounds  on  the  error  estimator,  after  the  last  iteration  step. 
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Figure  4.7:  Bandlimited  wave,  u>max  =  1000  Hz:  Comparison  of  the  exact  solution  with  the 
homogenized  solution  (left),  and  the  adapted  solution  (right);  all  are  normalized  by 
maximum  ofEdu(x,  t ) /dx,  x  €  fl,  t  G  [0, 0.057]. 
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In  Figures  4.8  through  4.10,  the  computed  evolution  of  stresses  is  shown.  As  ex¬ 
pected,  the  effect  of  dispersion  is  observed  to  increase  as  <umax  increases.  As  a  con¬ 
sequence,  the  homogenized  solutions  represent  a  very  poor  approximation  as  wmax 
reaches  higher  values.  The  phenomenon  of  trailing  waves  [3]  is  observed  in  particular 
for  the  solutions  for  camax  =  4000  Hz.  The  solutions  obtained  through  the  adaptive 
modeling  process  not  only  follow  the  exact  solution  in  the  domain  of  interest,  but  also 
globally  are  much  more  accurate. 
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Figure  4.8:  Bandlimited  wave,  u>max  =  2000  Hz:  Comparison  of  the  exact  solution  with  the 
homogenized  solution  (left),  and  the  adapted  solution  (right);  all  are  normalized  by 
maximum  of  Edu(x,t)/dx,  x  €  fl,  t  G  [0, 0.057]. 
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Figure  4.9:  Bandlimited  wave,  u>max  =  3000  Hz:  Comparison  of  the  exact  solution  with  the 
homogenized  solution  (left),  and  the  adapted  solution  (right);  all  are  normalized  by 
maximum  of  Edu(x,t)/dx,  x  €  fl,  t  G  [0, 0.057]. 
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Figure  4.10:  Bandlimited  wave,  wmax  =  4000  Hz:  Comparison  of  the  exact  solution  with  the 
homogenized  solution  (left),  and  the  adapted  solution  (right);  all  are  normalized  by 
maximum  of  Edu(x,t)/dx,  x  £Cl,t  £  [0,0.057]. 
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5  Concluding  Remarks 


A  mathematical  theory  and  computational  procedures  are  presented  for  a  systematic 
analysis  of  wave  propagation  in  elastic  heterogeneous  media.  The  issue  of  mathemat¬ 
ical  modeling  of  this  phenomenon  is  addressed  by  applying  a  general  notion  [17]  of 
residual-based  a  posteriori  analysis  of  errors  occurring  in  problems  in  computational 
engineering  sciences.  The  main  conclusions  of  this  work  can  be  summarized  as  fol¬ 
lows: 


•  The  general  framework  [17]  for  error  analysis  is  applied  to  derive  an  a  posteriori 
estimate  of  the  modeling  error  in  a  quantity  of  interest,  the  average  stress  on  a 
subdomain.  This  makes  possible  the  accurate  assessment  of  the  error  caused  by 
the  occurrence  of  geometric  wave  dispersion,  characteristic  to  dynamic  behav¬ 
ior  in  heterogeneous  media,  but  also  of  errors  caused  by  local  features,  such  as 
inaccurate  descriptions  of  the  local  material  model. 

•  Lower  bounds  to  a  global  norm  of  the  modeling  error  are  derived,  that  in  the 
presence  of  sufficient  damping  in  the  system,  represent  reasonably  accurate  a 
posteriori  estimates. 

•  Nonlocal  error  indicators  are  introduced  that  enable  nonlocal  model  adaptation 
and,  therefore,  make  possible  the  successful  control  of  geometric  dispersion. 

•  The  concept  of  goal-oriented  adaptive  modeling,  introduced  by  Zohdi,  Oden 
and  Rodin  [29],  and  Oden  and  Vemaganti  [19,  20,  26,  27],  is  extended  and  suc¬ 
cessfully  applied  to  the  elastodynamic  wave  problem.  An  adaptive  modeling 
algorithm  is  presented  that  solves  the  wave  problem  in  the  frequency  domain 
and  allows  an  adaptation  of  the  models  only  for  those  waves  that  contribute 
most  to  the  error.  In  practical  applications,  this  implies  that  the  algorithm  makes 
possible  the  identification  of  waves  that  are  most  dispersion  sensitive  to  the  di¬ 
mension  of  the  inhomogeneity  in  the  material.  The  error  indicators  can  then  be 
used  to  adapt  the  material  model  nonlocally  for  these  critical  frequencies  and 
enable  the  identification  of  those  regions  in  the  elastic  body  that  contribute  to 
dispersion. 

As  future  extensions  of  this  work,  one  could  consider  the  application  of  the  adap¬ 
tive  modeling  methodology  to  the  inverse  problem,  i.e.  knowing  the  response  one 
would  seek  the  actual  material  model  within  a  degree  of  certainty  or  accuracy.  This 
application  finds  its  interest  in  many  disciplines  of  engineering  sciences,  such  as  e.g. 
seismology  and  structural  acoustics. 


67 


Also,  the  incorporation  of  control  of  the  numerical  approximation  error  into  the 
algorithm  is  a  possible  extension.  The  resulting  methodology  would  consequently 
address  both  the  issues  of  Verification  and  Validation  of  wave  propagation  in  hetero¬ 
geneous  media. 

Apart  from  these  areas,  there  are  many  fields  to  be  explored  for  future  extension 
of  the  GOALS  methodology,  such  as,  for  example:  nonlinear  elasticity,  visco-elasticity, 
plasticity,  and  fluid  flow  problems. 
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Appendix  A 

Derivation  of  Optimal  Lower  Bounds  //low( ;) 


In  this  section,  a  more  detailed  derivation  is  presented  of  the  bounds  77 jj^  that  are 
introduced  in  Lemma  3.4.1.  Since  the  derivations  of  the  four  bounds  are  similar,  the 
derivation  of  only  one  bound,  77^  ,  is  shown.  First,  by  recalling  the  result  in  the  proof 
of  Lemma  3.4.1,  one  gets: 


ei  +  eillw  > 


|7e(uo  +  po,uo  +  (9+po)| 

Huo  +  tf+pollw 


\/6+  eC. 


^low  chosen  to  be  equal  to  the  RHS  with  a  value  of  9  +  for  which  this  RHS  reaches  a 
maximum.  Clearly  a  lower  bound  that  is  closer  to  ||ei  +  £i||w,  provides  a  more  accu¬ 
rate  estimate  of  this  norm  and,  inherently,  the  modeling  error.  Hence,  the  RHS  has  to 
be  minimized  with  respect  to  9+ .  By  taking  the  following  expansion:  0+  =  a  +  hi,  the 
RHS  of  the  above  expression  can  be  rewritten  as  follows: 


|fc(u0  +po,u0  +  (a  +  6i)po)| 

||uo  +  (a  +  6i)p0||w 


<F(a,6), 


Va,6  G  M. 


By  applying  the  definition  (3.13)  of  the  norm  ||  •  \\ji  and  the  norm  definition  of  a  complex¬ 
valued  number:  \z\  =  \f//,  one  gets: 


^  / ^(u0  +  p0,u0+(a  +  bi)p0)K(u0  +  p0,uo+(a  +  bi) p0)  /ttttv 

®{a,  b)  \  ^  ,  1  i  ~  .  /  |  1  -  \  *  \  v  b) ' 

V  H(uo  +  (a  +  bi)po,U-o  +  (a  +  bi)Vo) 


The  minimizers  of  this  terms  are  numbers  a  and  b  such  that 
plies: 


<9$ 

da 


=  0.  This  im- 

db 


<9<F  d(p 
dp  da 

<9<F  dp 
dp  db 


1  dp 

‘l  s/p  da 

1  dp 
2  s/p  db 


0, 

0, 


dp 

da 


=  0, 


dp 

~db 


=  0. 
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dp 

After  several  algebraic  manipulations,  the  following  is  obtained  for  — : 


dp 

da 


|u0  +  (a  +  6i)p0||^ 


He  |7e(uo  +  po,  Po)^(«o  +  Po,  uo) }  |]uo \\] 
-|^(u0  +  p0,uo)|2  He{H{ u0,po)} 


+2  a 


|7e(u0  +  po,po)|2  ]|u0||w  -  |^(u0  +  p0,uo)|2  llpol 


n 


+4  ab 


I H 


Tm  {7e(uo  +  po,  uo)^(uo  +  Po,  Po) }  ||po 

+|7e(u0  +  po,  Po)  |2  Irn {H{ u0,  p0)} 


—2a2 


I  n 


+46 


Tie  |7?-(u0  +  p0,uo)7\!.(u0  +  po,po)|  ||po 

-|7e(u0  +  po,po)|2  He {7t(u0,p0)} 

He  {ft(u0  +  po, po)^(u0  +  po, uo) }  lm {7t(u0, po)} 

+lm  [h(vl0  +  po, u0)^-(uo  +  Po,  Po) }  He  {H(u0,  po)} 

+262  |  He  {ft(u0  +  po,  po)^(uo  +  po,  uo) }  ||po \\n 

-|^(u0  +  p0,po)|2  He{H(u0,p0)} 
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dp 

Whereas,  for  — -  one  gets: 
ob 


dp 

~db 


|u0  +  (a  +  6i)po||^ 


1m  {ft(u0  +  p0,  po)^(uo  +  Po,  u0) }  ||u0 \\] 
-|^(u0  +  po,u0)|2  Jm{H(u0,po)} 


+4  a 


1m 


{n(u0  +  Po, po)^(u0  +  Po, Uo) }  Tie {H(u0, Po)} 

—Tie  {ft(u0  +  po,  u0)^(u0  +  po,  po) }  Trn {H(u0,  p0)} 


—4  ab 


I H 


^e|^(u0  +  po,u0)^(u0  +  po,Po)}  ||po 

-|^(u0  +  p0,Po)|2  ^e{7f(u0,p0)} 


+2  az 


I  n 


J?n|7e(uo  +  po,Po)^'(uo  +  po,uo)|  ||po 

-|7e(u0  +  po, Po) |2  lm{H(u o,po)} 


+26 


+262 


|^(u0+p0,po)|2  ||u0||w  -  |^(u0  +  po,u0)|2  ||po| 


n 


I  n 


lm  |7^.(u0  +  po,  u0)^.(u0  +  po,  po)  |  ||po 

+  |^(u0  +  p0,po)|2  lm{H{ u0,po)} 


To  simplify  notations,  pa  and  pi,  are  introduced,  such  that: 

dp  _  1  „  dp  _  1 

da  ||u0  +  (a  +  6?')p0||^  ^a'  db  ||u0  +  (a  +  6i)p0||^  ^b' 

It  is  clear  that  a  and  6  that  maximize  the  lower  bound,  satisfy  pa  =  Pb  =  0.  This  non¬ 
linear  set  of  equations  is  solved  numerically  by  applying  a  Newton-Raphson  scheme. 
Thus,  an  iterative  procedure  is  started,  where  at  step  i  the  values  for  a  and  6  are  de¬ 
termined  by  ai  =  aj_i  +  A  a,  and  6,  =  6*_i  +  A6„  where  A  a,  and  A  6,  are  computed  by 
solving: 


/  (  dpa\ 

V  da 

(  dTh>\ 

\\da  yi_1 


dpa\  \ 

db  A-r 

dp  b\ 

db  A-i/ 


^  A  a, 

yAbi 


\(<Pb)i- 1, 
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where  the  coefficients  of  the  matrix  are  given  by: 

~^=  2  l^(uo  +  Po,Po)|2  lluollw  -  |^(uo  +  Po,u0)|2  Upollw 
+46  Tm|^(u0  +  po,Uo)^(uo  +  po,po)}  Hpollw 

+|ft(u0  +  po,po)|2Xra{?f(uo,po)} 
-4 a  7ee|7e(uo  +  po,u0)^(u0  +  po,po)}  llpollw 

-|7e(u0  +  po,Po)|2^e{?f(uo,po)}  , 

-^=  4a  J?n|^(u0  +  po,uo)^(u0  +  po,po)}  IIpoIIk 

+|^(u0  +  p0,po)|22:m{W(uo,po)} 

+4  ^e|7e(u0  +  po,po)^(uo  +  Po,uo)}  {7+(u0, Po)} 
+lm  |^(u0  +  po,  uo)^(uo  +  Po,  Po) }  Tie  {H{ u0,  po)} 

+46  7ee|^(u0  +  po,po)^(uo  +  Po,uo)}  Hpollw 

-|^(u0  +  p0,po)|2  Tie{H(u0,po)}  , 

~j^  =  4  2Tto|^(u0  +  Po,Po)^(uo  +  Po,Uo)}  TZe{H{ u0,Po)} 
-lZe^lZ(ixQ  +po,u0)^(u0  +p0,po)}  Xm{7f(u0,p0)} 

-46  ^e|^(uo  +  po,u0)^(u0  +  po,po)}  llpollw 

-|7e(u0  +  po,Po)|2  Tie{H{ u0,Po)} 
+4a  Xm|^(u0  +  po,po)^(uo  +  po,u0)|:  Hpollw 

—  |^(u0  +  POiPo)|2  2Tm{7i(uo,po)}  , 
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~^  =  ~4a  ^e|^(u0  +  po,uo)^(u0  +  Po,Po)}  llpollw 

-|^(u0  +  po,po)|2  ^e{H{u0,po)} 

+2  |7e(u0+p0,po)|2  IIuoIIk  -  |^(uo  +  Po,uo)|2  llpollw 
+46  T?n|^(u0  +  po,uo)^(uo  +  Po,Po)}  ||po||« 

+|7e(u0  +  p0,po)|22:^{W(uo,Po)}  ■ 
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